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TIME  SERIES  DATA  ANALYSIS 
AND  SYNTHESIS  FOR  RESEARCH  WATERSHEDS 


By  W.  M.  Snyder1 

ABSTRACT 

Mathematical  models  of  watershed  flow  processes  can  be  used  to  simulate 
future  runoff  events.  Such  synthesized  series  of  events  form  the  basis  for  risk 
analysis  in  planning  and  design.  Operation  of  the  models  requires  generation  of 
synthetic  rainfall  and  other  model  inputs  as  stochastic  series.  This  paper  presents 
model  forms  for  analysis  and  synthesis  of  information  usually  available  from 
research  watersheds. 


INTRODUCTION 

Fiering  and  Jackson2  have  presented  the  practi- 
cal approach  to  generation  of  synthetic  streamflows 
in  an  excellent  and  understandable  tutorial  mono- 
graph. A  series  of  five  computer  programs  has  been 
designed  to  implement  such  generating  schemes. 
Four  of  the  programs  are  designed  to  examine  his- 
torical data  as  an  aid  in  choosing  from  among  simple 
linear  generating  models.  The  fifth  synthesizes 
streamflows  for  a  particular  model. 

The  programs  are  presented  in  terms  of  case 
studies  on  two  watersheds  operated  by  the  Agricul- 
tural Research  Service — Watkinsville,  Ga.,  W-l, 
19.2  acres;  and  Vero  Beach  (Taylor  Creek),  Fla., 
W-2,  98.6  square  miles.  Data  for  these  watersheds 
are  published  in  Hobbs3  and  subsequent  reports. 
Published  data  for  1940-64,  were  available  for  the 
Watkinsville  watershed.  Data  for  1956-66  were 
available  for  Taylor  Creek  watershed.  In  addition, 
mean  monthly  air  temperatures  at  Lake  Okee- 
chobee, Fla.,  and  Athens,  Ga.,  were  taken  from 
Climatic  Data  publications  of  the  U.S.  Wea- 
ther Bureau  (now  Weather  Service).  No  fur- 


1Hydraulic  engineer,  Southeast  Watershed  Laboratory,  Ag- 
ricultural Research  Service,  U.S.  Department  of  Agriculture, 
Athens,  Ga.  30601. 

2Fiering,  M.  B.,  and  Jackson,  B.  B.  1971.  Synthetic  stream- 
flows.  Am.  Geophys.  Union  Water  Resour.  Monogr.  1,  98  pp. 

3Hobbs,  H.  W.  1963.  Hydrologic  data  for  experimental  ag- 
ricultural watersheds  in  the  United  States.  U.S.  Dep.  Agric. 
Misc.  Publ.  No.  945,  pp.  8.2-1-8.2-4,  10.1-1-10.1-8. 


ther  descriptions  of  the  watersheds  or  the  data  will 
be  given,  since  this  report  deals  with  development 
of  numerical  methodologies,  not  with  analysis  of  any 
particular  data  set. 

Data  from  the  two  Agricultural  Research  Service 
watersheds  were  used  so  that  implementation  of  the 
methods  and  models  suggested  by  Fiering  and 
Jackson  could  be  tested  on  watersheds  smaller  than 
those  normally  considered  for  water  resource  de- 
velopment. Such  testing  is  a  necessary  part  of  the 
development  of  procedures  to  generate  synthetic 
inputs  for  conceptual  watershed  models.  Current 
emphasis  in  agricultural  watershed  modeling  is  on 
prediction  of  water  runoff  as  the  earner  of  sediment 
and  other  materials.  Ranges  and  extremes  of  water, 
sediment,  and  chemical  runoff  must  be  estimated  by 
techniques  of  synthetic  data  generation. 

This  report  presents  some  results  of  application 
of  numerical  methods  not  covered  by  Fiering  and 
Jackson.  Regression  coefficients  are  estimated  by 
both  conventional  and  nonlinear  least  squares. 
Statistical  frequency  distributions  with  mathemati- 
cally continuous,  seasonally  variable  parameters 
are  also  estimated  by  nonlinear  least  squares.4  Sea- 
sonally continuous  regression  coefficients  and  dis- 


"Snyder,  W.  M.  1972.  Fitting  of  distribution  functions  by 
nonlinear  least  squares.  Water  Resour.  Res.  8(6):  1423-1432, 
Snyder,  W.  M.,  and  Wallace,  J.  R.  1974.  Estimating  the 
parameters  of  the  log-normal  distribution.  Nordic  Hydrol. 
5(3):  129-145. 
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tribution  parameters  are  specified  by  a  cyclic  mod- 
ification of  the  method  of  continuous  parabolic 
interpolation.5 

In  the  following  presentation,  no  particular  at- 
tempt is  made  to  interpret  the  hydrologic  implica- 
tions of  all  outputs  from  the  computer  programs, 
but  rather  to  illustrate  how  the  programs  provide  a 
basis  for  decisionmaking  on  model  structure,  as  dis- 
cussed by  Fiering  and  Jackson.  The  programs  are 
not  intended  solely  for  use  in  their  listed  form.  They 
do  form  core  procedures  which  can  be  easily  mod- 
ified to  incorporate  other  operations  on  data  and 
other  outputs  as  required  by  individual  research- 
ers. For  example,  many  models  will  be  based  on 
other  than  monthly  data,  and  modifications  can  be 
made  to  the  programs  to  change  the  time  base  of 
the  data. 

For  brevity  of  presentation,  it  is  assumed  that  the 
reader  is  familiar  with  time-series  analysis  and 
synthesis  from  tutorial  and  applications  perspec- 
tives such  as  presented  by  Fiering  andJackson. 


AUTOCORRELATION  AND 
CROSS  CORRELATION 

Consider  the  two  series  of  variables:  Yh  i=l,  N; 
andXh  i  =  l,N.  Here,  Y,  through  YN  will  represent 
monthly  totals  (or  averages)  of  some  quantity,  such 
as  watershed  runoff,  in  chronological  order  from 
month  1  to  month  N  of  a  record.  Similarly,  X/s  are 
such  totals  or  averages.  The  Z,'s  represent  inde- 
pendent, or  causal,  properties.  The  F(-'s  represent 
dependent  properties. 

The  end  product  of  a  data  generation  scheme  is 
the  estimation  of  one  or  more  possible  sets  of 
the  series  Y{.  For  example,  one  50-year  series  of 
monthly  runoff  values  might  be  synthesized,  or  10 
such  series  might  be  synthesized.  The  synthetic 
series  are  estimates  of  what  might  happen  in  the 
future.  Assurance  is  needed  that  these  predicted 
series  are  likely,  to  some  degree  of  likeliness,  usu- 
ally based  on  one  observed  series  of  Yt. 

The  likeliness  of  possible  future  series  of  values  of 
Y,  is  accomplished  by  extracting  certain  elements  of 
information  from  the  historical  sequence  and  mak- 
ing certain  that  this  same  information  is  reproduced 
in  the  synthetic  series.  For  example,  the  means  and 
standard  deviations  of  many  synthetic  series  should 
have  mean  values  very  nearly  the  same  as  the  his- 
torical series. 

'Snyder,  W.  M.  1961.  Continuous  parabolic  interpolation. 
J.  Hydraul.  Div.,  Proc.  Am.  Soc.  Civ.  Eng.  87(HY4):  99-111. 


The  starting  point  for  a  data-generating  scheme 
must  be  a  method  of  processing  data  so  that  infor- 
mation within  the  data  can  be  recognized  and  quan- 
tified. It  must  be  determined  whether  the  Y-s  are 
serially  related.  Is  there,  for  example,  a  relation- 
ship between  Yi=p  and  Yi=p+l,  where  I  is  some  posi- 
tive number  of  lags?  Does  this  relation  vary,  and  in 
particular  does  it  tend  to  have  a  cyclic  pattern  as  / 
increases  from  1  to  some  number  less  than  N1  Such 
cyclic  within-series  patterns  of  information  are  usu- 
ally determined  by  autocorrelation,  which  simply 
means  the  calculation  of  the  correlation  coefficient 
between  Yi=p  and  Yi=p+l  for  various  values  of  I. 

When  a  supplemental  data  series,  X, ,  is  available, 
relationships  between  this  series  and  the  series  to 
be  predicted  should  also  be  determined.  If  monthly 
rainfall  is  available,  synchronous  with  the  monthly 
runoff  values,  then  it  must  be  determined  whether 
information  in  the  runoff  series  can  be  cast  into  a 
rainfall-runoff  relationship.  Also,  it  must  be  deter- 
mined whether  this  relationship  is  more  important 
than  the  serial  structure  of  the  Y,  series  in  the 
prediction  of  F,  . 

Program  I,  listed  in  the  appendix,  is  a  simple 
program  designed  to  generate  the  autocorrelation 
coefficients  for  one  independent  and  one  dependent 
variable  and  for  the  cross-correlation  coefficients 
between  the  two.  The  dependent  variable  lags  be- 
hind the  independent  variable,  since  it  is  obvious, 
for  example,  that  runoff  cannot  anticipate  future 
rainfall. 

Figure  1  is  a  plot  of  output  from  program  I  for  the 
Taylor  Creek  watershed.  In  the  upper  part  of  the 
figure,  the  auto-  and  cross-correlations  for  monthly 
rainfall  and  runoff  are  shown.  In  the  lower  part,  the 
autocorrelation  coefficient  for  monthly  mean  tem- 
perature and  cross-correlation  coefficients  between 
temperature  and  rainfall  and  runoff  are  shown. 
Temperature  is  considered  the  independent  vari- 
able. The  coefficients  were  computed  for  all  lags 
from  0  to  39. 

The  strong  seasonal  pattern  of  rainfall,  runoff, 
and  temperature  is  shown  by  the  sinusoidal  waves 
with  a  wave  length  of  12  lags  in  both  the  upper  and 
lower  portions  of  figure  1.  Temperature  shows  the 
most  distinct  pattern,  with  correlations  oscillating 
between  values  of  -0.9  and  +0.9.  Rainfall  also 
shows  a  distinct  annual  pattern,  but  values  oscillate 
between  about  -0.4  and  +0.4,  showing  a  much 
greater  year-to-year  variability  than  temperature. 
Monthly  runoff  has  the  lowest  autocorrelation  co- 
efficients, and  hence  the  greatest  year-to-year 
variability. 
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Figi  rk  1. — Auto-  and  cross-correlation  coefficients  for  Tavlor  Creek  watershed. 
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The  seasonal  patterns  of  the  cross-correlation 
coefficients  show  that  rainfall,  runoff,  and  tempera- 
ture are  all  in  phase,  highs  of  rainfall  and  runoff 
occurring  with  highs  of  temperature.  This  is  an 
expression  of  the  semitropical  climate  of  Florida, 
where  the  summer  season  is  the  wet  season  of  the 
year. 

Figure  2  is  a  plot  of  the  output  from  program  I  for 
the  Watkinsville  watershed.  In  this  figure  it  is  read- 
ily apparent  that  temperature  follows  the  same 
strong  seasonal  pattern  as  shown  in  Florida.  How- 
ever, rainfall  and  runoff  show  only  indistinct  sea- 
sonal patterns.  The  lower  part  of  figure  2  shows  that 
rainfall  and  runoff  for  this  watershed  are  not  in 
phase  with  temperature.  In  a  temperate  climate 
under  variable  continental  and  maritime  influence, 
the  wet  season  tends  to  be  in  late  winter  and  early 
spring. 

Figures  1  and  2  indicate  that  both  rainfall  and 
runoff  have  higher  serial  correlations  with  tempera- 
ture than  they  do  with  themselves  or  with  each 
other.  It  is  evident  that  a  procedure  for  synthesiz- 
ing a  runoff  series  must  take  into  account  this  strong 
seasonal  pattern.  Temperature  could  be  used  as  an 
index  of  the  seasonal  pattern.  However,  since  the 
autocorrelation  coefficients  for  temperature  show 
a  seasonally  repeating  pattern  of  about  0.9,  the 
month  of  the  year  is  also  a  good  indication  of  season 
and  is  much  easier  to  use. 


AUTOREGRESSION  AND 
CROSS-REGRESSION 

The  information  in  figures  1  and  2  implies  a  re- 
gression of  the  type 

ROxl=a+b1Pu+c1ROu-1+b,Pu-,+c<,ROu-,  +  . . . 

(1) 

In  equation  1,  ROM  is  runoff  for  month  M,  PM  is 
precipitation  for  month  M,  and  ROM^  and  PM-\  are 
runoff  and  precipitation  for  the  antecedent  month. 
Additional  antecedent  months  can  be  included.  The 
terms  a,  b,  and  c  are  regression  coefficients.  Rain- 
fall and  runoff  during  any  month  are  significantly 
correlated,  and  therefore  it  should  not  be  necessary 
to  include  both  back  rainfall  and  back  runoff  as 
functional  (not  statistical)  independent  variables. 

Equation  1  would  usually  be  simplified  to 

ROM=a+b1PM+C1ROM-1+C2ROM_2+C3ROM-3+ . . 

(2) 

Equation  2  implies  an  operating  position  for  pre- 
diction as  follows:  A  rainfall  record,  Pu,  is  available 


by  months,  and  values  of  runoff  for  some  number  of 
months  prior  to  the  first  month  of  rainfall  record  are 
also  available.  Runoff  can  then  be  predicted  for  the 
first  month  of  the  rainfall  record.  Using  this  syn- 
thesized runoff  as  back  runoff  with  the  rainfall  for 
the  second  month,  runoff  for  the  second  month  can 
be  generated;  the  process  is  then  repeated  for  the 
full  length  of  a  real  or  synthetic  rainfall  record,  PM. 

Small  watersheds  may  be  defined  as  those  drain- 
age areas  for  which  soil-water  processes  are  pre- 
dominant over  channel  processes.  For  such  water- 
sheds the  number  of  back  values  of  runoff  that  can 
influence  runoff  of  a  succeeding  month  must  be 
small.  The  number  of  back  months  can  be  estimated 
from  graphs  such  as  figure  1  or  2.  The  number  is  the 
number  of  lags  before  auto-  or  cross-correlation 
first  drops  to  near  zero.  Fiering  and  Jackson6  (p.  67) 
give  a  systematic  test,  based  on  the  coefficient  of 
determination,  to  aid  in  choosing  the  number  of 
lags.  Program  II  in  the  appendix  uses  only  one  back 
month,  as  in 

ROM=a+bPM+cROM^.  (3) 

It  would  be  a  simple  matter  to  modify  this  pro- 
gram to  include  additional  months  by  using  a  ma- 
trix-inversion library  routine. 

It  is  usually  necessary  to  let  the  final  generating 
method  evolve  from  systematic  examination  of  the 
data  set.  Therefore,  program  II  is  constructed  to 
evaluate  all  the  regression  coefficients  in  the  set  of 


six  equations,  4  through  9. 

ROv=a+cROM_x  (4) 

RO„=a+bP„  (5) 

ROu=a+bP„+cROM-x  (6) 

\ogROu=a  +c  log  RO,,^  (7) 

log ROM=a+  b  log  Pu  (8) 

\ogROxl=a+b  logPu+c  log  RO,,_l  (9) 


Each  of  these  equations  is  evaluated  separately 
for  each  calendar  month  of  the  year;  the  choice  of 
model  can  thus  be  guided  by  best-fit  considerations 
as  well  as  operational  requirements. 

Tables  1  through  4  show  the  coefficients  derived 
by  least-squares  fitting  of  the  equations  to  the  data 
sets  for  Taylor  Creek  and  Watkinsville.  Cases  A,  B, 
and  C  for  natural  data  refer  respectively  to  equa- 
tions 4,  5,  and  6.  Cases  A,  B,  and  C  with  trans- 
formed data  refer  to  equations  7,  8,  and  9. 

Figure  3  is  a  plot  of  the  coefficients  of  equation  6 
for  Taylor  Creek.  Figure  4  is  a  similar  plot  for  the 
Watkinsville  watershed.  The  lower  part  of  each  of 
these  figures  shows  the  residual  errors,  consisting 


6Cited  in  footnote  2. 
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Figi'RK  3.— Monthly  regression  coefficients  and  residual  errors  for  Taylor  Creek  watershed. 
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Figure  4.— Monthly  regression  coefficients  and  residual  errors  for  Watkinsville  watershed. 
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TABLE  1. — Taylor  Creek  natural  data:  Auto-  and  cross-regression  coefficients  by  month 


Case  A1  Case  B2  Case  C3 


iviunin 

-  . ... 

a 

c 

a 

& 

a 

b 

c 

A  OOQ 

0.935 

a  i  a~ 
—  U.  iUo 

0.270 

A  OAT 

—  0.601 

A  OTC 

1.101 

A  Q1 

1  AA 

—  .Zoo 

0 1  o 

A  AC 

—  .44b 

OO  1 

.oo4 

OTA 

QQ7 

1  rtOO 

—  .bZo 

A  TO 
.475 

OAyl 

—  .o94 

1  oo 

.819 

1  A  Q 

/17Q 

.U7o 

AO  A 

—  .064 

1  AO 

.10Z 

1  oo 

—  .loo 

1  AA 

ACQ 

Mav   

AC1 

. .  .  —.Obi 

i  ai  i 

—  .OlZ 

Ol  A 
.^19 

™1  o 
—  .016 

.  1^4 

l.obl 

June   

i  cne 

. . .  l.byb 

—  .114 

—2.771 

.579 

—3.308 

.616 

.710 

T.,1,. 

1  171 

1  OA 

—  Z.411 

.DO  1 

O    1  AT 

— ZAOI 

.boo 

AAC 

.Uvb 

937 

.  U  I  u 

—  1  9^7 

.  OU  I 

.  t70W 

.oo*. 

.422 

Sept  

1.000 

1.067 

-2.365 

.802 

-2.598 

.676 

.535 

Oct  

1.531 

.195 

-.636 

.699 

-2.401 

.805 

.418 

Nov  

.177 

.061 

-.060 

.289 

-.245 

.310 

.073 

Dec  

...  .116 

.149 

.141 

.014 

.090 

.167 

.151 

Equation  4. 
2Equation  5. 
3Equation  6. 


TABLE  2. — Taylor  Creek  log -transformed  data:  Auto-  and  cross-regression  coefficients  by  month 
„     ,  Case  A1  CaseB2  CaseC3 


iviontn 

a 

c 

a 

b 

a 

b 

c 

Jan  

-0.359 

0.547 

-1.935 

0.999 

-0.883 

0.974 

0.495 

Feb  

-.248 

.815 

-2.364 

1.168 

-1.067 

1.015 

.728 

Mar  

-.089 

.811 

-2.831 

1.536 

-1.547 

1.561 

.824 

Apr  

-1.337 

.600 

-2.728 

.702 

-1.932 

.747 

.605 

May   

-1.005 

.410 

-3.465 

1.156 

—2.575 

1.145 

.405 

June   

.263 

.469 

-7.190 

3.335 

-6.281 

3.302 

.446 

July   

.161 

.643 

-7.231 

3.888 

-5.646 

3.101 

.271 

Aug  

-.009 

1.029 

-5.875 

3.075 

-3.221 

1.732 

.837 

Sept  

.531 

.437 

-3.559 

2.178 

-3.005 

1.917 

.302 

Oct  

-.440 

.678 

-1.069 

.832 

-1.528 

.969 

.734 

Nov  

-1.678 

.420 

-1.752 

.361 

-1.682 

.552 

.478 

Dec  

. . .  .-1.207 

.502 

-2.127 

.396 

-1.316 

.334 

.461 

'Equation  7. 
2Equation  8. 
3Equation  9. 


Table  3. 

— Watkinsville  natural  data:  Auto-  and  cross-regression  coefficients  by  month 

Month 

Case  A1 

Case  B2 

Case  C3 

a 

c 

a 

b 

a 

b 

c 

Jan. 

0.164 

1.410 

-1.021 

0.326 

-1.014 

0.272 

1.065 

Feb. 

.444 

-.061 

-.482 

.195 

-.552 

.203 

.060 

Mar. 

.417 

.522 

-1.440 

.349 

-1.451 

.337 

.193 

Apr. 

.243 

.457 

-.632 

.265 

-.695 

.233 

.327 

Mav 

.053 

.302 

-.202 

.118 

-.432 

.134 

.331 

June 

.272 

-.105 

-.774 

.281 

-.771 

.281 

-.011 

July  . 

.439 

.143 

-1.210 

.334 

-1.214 

.333 

.030 

.422 

-.077 

-.734 

.302 

-.735 

.303 

.002 

Sept. 

.012 

.033 

.005 

.007 

-.036 

.015 

.040 

Oct. 

.061 

.553 

-.102 

.069 

-.107 

.068 

.337 

Nov. 

.352 

-.513 

-.822 

.346 

-.921 

.358 

.792 

Dec.  . 

.  .254 

-.037 

-.406 

.147 

-.394 

.146 

-.031 

'Equation  4. 
2Equation  5. 
3Equation  6. 
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TABLE  4. — Watkinsville  log -tram formed  data:  Auto-  and  cross-regression  coefficients  by  month 


Case  A1  Case  B2  Case  C3 


a 

c 

a 

b 

a 

6 

c 

 -1.416 

0.374 

-7.532 

3.402 

-6.556 

3.258 

0.238 

Feb  

 -1.825 

.183 

-5.819 

2.478 

-5.306 

2.699 

.320 

Mar  

 -1.209 

.318 

-8.656 

3.922 

-7.914 

3.863 

.279 

Apr  

 -1.763 

.242 

-4.768 

1.965 

-4.337 

1.935 

.202 

May  

 -3.228 

.130 

-5.057 

1.364 

-4.765 

1.365 

.131 

June   

 -2.785 

.219 

-5.894 

2.027 

-5.163 

2.014 

.203 

July   

 -.654 

.669 

-6.673 

2.437 

-4.248 

2.004 

.500 

Aug  

 -2.902 

.180 

-5.411 

1.822 

-4.952 

1.802 

.144 

Sept  

 -3.439 

.187 

-4.595 

.532 

-3.884 

1.006 

.338 

Oct  

 -3.115 

.225 

-4.333 

.593 

-3.475 

.589 

.210 

Nov  

 -4.143 

-.097 

-5.150 

1.515 

-3.957 

1.815 

.364 

Dec  

 -3.212 

-.020 

-7.180 

2.911 

-7.604 

2.938 

-.103 

Equation  7. 
2Equation  8. 
3Equation  9. 


of  the  difference  between  the  observed  runoff  and 
the  predicted  runoff  for  each  month  of  record;  they 
are  grouped  by  classes.  Predicted  runoff  is  gener- 
ated with  the  coefficients  to  appropriate  to  the 
calendar  month,  as  given  in  the  upper  part  of  the 
figure. 

The  coefficients  in  figures  3  and  4  exhibit  in  a 
crude  way  the  seasonal  patterns  of  figures  1  and  2. 
The  pattern  is  irregular  because  the  statistical  av- 
eraging in  figures  3  and  4  is  over  only  approximately 
one-twelfth  of  the  data  points  in  figures  1  and  2.  For 
example,  a  lag-1  correlation  coefficient  is  based  on 
all  pairs  January-February,  February-March, 
March-April,  and  so  on  through  the  years  of  record. 
The  monthly  regression  coefficients  for  a  calendar 
month,  however,  are  based  only  on  observations  for 
that  calendar  month  in  the  record. 

SEASONALLY  CONTINUOUS 
AUTOREGRESSION  AND 
CROSS  REGRESSION 

This  section  describes  the  development  of  a  mod- 
ified form  of  equation  6  that  is  mathematically  con- 
tinuous in  a  cyclic  pattern  throughout  the  year.  It 
should  be  remembered  that  the  regression  coeffi- 
cients in  figures  3  and  4  are  in  discrete  sets,  each 
calendar  month  producing  a  set.  A  total  of  36  coeffi- 
cients are  needed,  and  must  be  obtained  from  the 
data  by  least  squares.  There  is  no  continuity  from 
month  to  month,  and  the  patterns  of  the  coefficients 
are  quite  irregular. 

Equation  6  can  be  modified  by  making  coefficients 
a,  b,  and  c  functions  of  calendar  month.  In  such  a 
modified  equation,  one  would  normally  be  required 


to  estimate  the  parameters  of  the  three  calendar- 
month  functions.  A  simpler  and  more  direct  method 
of  establishing  continuity  of  the  coefficients  is  by 
fitting  an  interpolation  function  by  least  squares. 
The  interpolating  function  used  here  was  continu- 
ous parabolic  interpolation7  modified  to  a  cyclic 
form. 

Consider  the  arrangement  of  monthly  values  of 
the  coefficient  a  of  equation  6,  shown  in  table  5. 
Base  values  are  located  as  shown,  one  every  4 
months.  Continuous  parabolic  interpolation  re- 


7Snyder,  cited  in  footnote  5. 

Table  5. — Schematic  for  seasonally  continuous 
interpolation 


Month 


Base 
value 


Interpolated 
value 


Sept. 

Oct.  . 

Nov. 

Dec. 

Jan.  . 

Feb. 

Mar. 

Apr. 

May  . 

June 

July  . 

Aug. 

Sept. 

Oct.  . 

Nov. 

Dec. 

Jan. 

Feb. 

Mai-. 

Apr. 

May 
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quires  four  base  values  for  an  interpolated  value. 
The  values  for  a2,  a3,  and a4  are  based  on  a9,  a,,  a5, 
and  a9.  It  can  be  seen  that  a  seasonally  continuous 
cyclic  system  of  interpolation  is  produced. 

Program  III,  listed  in  the  appendix,  was  written 
to  derive,  simultaneously  by  least  squares,  the  base 
values  a  u  a5,  a9,bu  bh,  b9,cu  c5,  andc9.  The  squares 
of  errors  for  every  month  of  the  record  are  utilized, 
rather  than  just  the  errors  for  those  months  where 
the  base  values  are  located.  It  can  be  noted  in  the 
program  listing  that  interpolation  is  earned  out 
using  numerical  interpolating  coefficients.  These 
coefficients  were  taken  from  table  7  in  Snyder.8 

Tables  6  and  7  are  listings  of  the  base  values, 
defined  in  table  5,  for  the  three  seasonal  functions 
for  the  coefficients  of  equation  6.  Several  different 
determinations  of  these  seasonal  functions  were 
made,  using  different  fitting  criteria.  The  criteria 
for  and  deviations  of  coefficients  are  discussed 
below. 

The  coefficients  for  natural  data,  unrestrained,  in 
tables  6  and  7  are  the  seasonally  continuous  versions 
of  the  discrete  sets  of  coefficients  in  figures  3  and  4. 
In  further  explanation,  in  table  8  the  coefficients  for 
natural  data,  unrestrained,  are  expanded  to  include 
all  calendar  months,  using  the  same  interpolating 
coefficients  as  in  the  listing  of  program  III.  These 
expanded  coefficients  are  plotted  as  continuous 

8Snyder,  W.  M.  1962.  Closure  to  continuous  parabolic  interpo- 
lation. J.  Hydraul.  Div..  Proc.  Am.  Soc.  Civ.  Eng.  88  (HY4): 
265-274. 


lines  in  figure  5  for  both  Taylor  Creek  and  Watkins- 
ville.  For  convenient  comparison,  the  discrete  coef- 
ficients from  figures  3  and  4  are  also  plotted.  It  can 
be  seen  that  the  cyclically  continuous  coefficients 
form  a  smooth  average  of  the  discretely  derived 
coefficients.  By  evaluating  only  9  parameters  from 
the  historical  record,  instead  of  the  36  in  monthly 
discrete  analysis,  a  greater  degree  of  averaging 
is  imposed.  This  averaging  tends  to  smooth  the 
irregularities  in  the  discrete  coefficients. 

The  comparisons  evident  in  figure  5  are  based  on 
equation  6  as  the  predictor  model.  Similar  compari- 
sons could  be  made  for  all  other  models  implied  in 
equations  4  through  9.  The  researcher  may  need  to 
make  such  additional  comparisons  before  choosing 
the  predictor  model  for  a  specific  watershed. 

Before  discussing  some  of  the  additional  fitting 
criteria  in  tables  6  and  7,  it  is  necessary  to  consider 
the  plots  in  figure  6.  The  upper  plot  shows  the 
position  a  line  would  take  if  fitted  by  customary 
least-squares  procedures.  In  particular,  the  posi- 
tive error  (observed  minus  calculated)  for  data  point 
1  would  be  roughly  balanced  by  the  negative  error 
for  point  3.  If  this  fitted  line  is  subsequently  used  for 
simulating  values  of  Y  based  on  values  ofZ,  a  bias 
results.  To  the  left  of  point  2,  calculated  values  of  Y 
are  negative.  These  can  be  set  to  zero,  since  nega- 
tive runoff  does  not  occur.  However,  from  point  2  to 
point  3,  calculated  values  of  runoff  are  positive, 
whereas  in  reality  they  should  also  be  zero.  Since 
the  positive  values  from  point  2  to  point  3  are  not 


TABLE  6. — Taylor  Creek  seasonally  cyclic  coefficients  in  ROM=aM+bMPM+cMROM-, 


Fitting  criteria 


Coeffi- 
cient 

C1) 

(2) 

Natural 
RO^O 

Natural 

6«1 
RO^O 

Estimated 
natural  for 
simulation 

a, 

-.6763 

-0.1540 

-12.1178 

-24.7986 

-3.50 

o5 

-3.0677 

-1.0969 

-3.6173 

-4.9686 

-1.71 

a„ 

-2.5040 

-2.1945 

-2.6206 

-3.2856 

-2.34 

b, 

.7140 

.3653 

2.2952 

1.0000 

.85 

I*: 

1.7859 

.3411 

.7240 

.8315 

.46 

b. 

1.5078 

.6773 

.6865 

.7178 

.68 

c, 

.6396 

.2367 

1.2124 

2.0082 

.49 

c3 

.5463 

.0302 

.1887 

.1050 

.07 

.5400 

.3246 

.4268 

.5433 

.35 

s„ 

.8819 

1.0965 

.9487 

1.5490 

1.05 

Y 

-1.1115 

1.1455 

1.0458 

1.1709 

1.1439 

Y 

-1.1115 

1.1439 

1.1439 

1.1439 

1  Log-transformation  of  data— unrestrained. 
2Natural  data — unrestrained. 


S,    Standard  deviation  of  residual  errors. 

Y_    Mean  of  calculated  values. 
Y    Mean  of  observed  values. 
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TABLE  7. — Watkinsville  seasonally  cyclic  coefficients  in  ROM=aM+bMPM+cMRO 


Coeffi- 
cient 

Fitting  criteria 

"Matilda  1 

(J) 

(2) 

RO>0 

a, 

-5.5750 

-0.8104 

-3.1206 

a.-. 

-5.0523 

-.8910 

-2.0072 

a,, 

-3.6056 

-.5464 

-5.0934 

b, 

2.5766 

.2556 

.5482 

b:, 

2.1564 

.2684 

.4032 

L 

0,, 

.7799 

.2355 

.7652 

c, 

.2144 

.1165 

.2612 

c5 

.2632 

.2769 

.4624 

C\< 

.2432 

-.0124 

.1258 

Sr 

1.2789 

.5868 

.4022 

Y 

-3.1350 

.3380 

.3104 

Y 

-3.1350 

.3380 

.3380 

'Log- transformation  of  data- 

-unrestrained. 

2Natural  data — unrestrained. 


S,.    Standard  deviation  of  residual  errors. 

Y  Mean  of  calculated  values. 

Y  Mean  of  observed  values. 


balanced  by  negative  values  from  point  1  to  point  2, 
the  average  of  the  simulated  values  is  higher  than 
the  average  of  the  observed  values. 


Table  8. — Exa  mples  of  complete  sets  of  seasonally 
cyclic  coefficients 


Watershed 

Coefficient 

and 

month 

a 

b 

c 

Taylor  Creek: 

Jan.  ...  

. . .  .-0.1540 

0.3653 

0.2367 

Feb  

....  -.0222 

.3306 

.1817 

Mar  

....  -.4293 

.3127 

.1095 

Apr  

....  -.7804 

.3151 

.0493 

May  

....-1.0969 

.3411 

.0302 

June   

....-1.4340 

.4150 

.0775 

July  

....-1.8322 

.5272 

.1700 

Aug  

....-2.1372 

.6304 

.2661 

Sept  

. . .  .-2.1945 

.6773 

.3246 

Oct  

. . .  .-1.8351 

.6381 

.3323 

Nov  

....-1.1839 

.5438 

.3120 

Dec  

....  -.5279 

.4383 

.2760 

Watkinsville: 

Jan  

. . .  .-0.8104 

0.2556 

0.1165 

Feb  

....  -.8534 

.2604 

.1675 

Mar  

....  -.8887 

.2653 

.2228 

Apr  

....  -.9050 

.2686 

.2676 

May  

....  -.8910 

.2684 

.2769 

June   

....  -.8204 

.2621 

.2018 

July  

....  -.7072 

.2515 

.1342 

Aug  

....  -.5997 

.2410 

.0410 

Sept  

  -.5464 

.2355 

-.0124 

Oct  

  -.5739 

.2370 

-.0103 

Nov  

  -.6518 

.2427 

.0239 

Dec  

....  -.7430 

.2499 

.0723 

A  modification  of  the  fitting  procedure  should 
significantly  reduce  this  bias.  In  the  lower  plot  in 
figure  6  is  shown  the  position  a  line  will  take  if 
negative  calculated  values  are  set  to  zero  during 
fitting  of  the  line.  Errors  at  points  1  and  2,  and 
perhaps  at  point  3,  become  zero.  The  line  takes  the 
position  shown,  minimizing  error  for  points  3,  4,  and 
5.  This  procedure  is  not  the  same  as  deleting  points 
with  zero  runoff  from  the  data  set  prior  to  fitting.  It 
is  not  known  ahead  of  time  whether  point  3  would 
have  a  negative  or  positive  calculated  value;  conse- 
quently, it  is  not  known  ahead  of  time  whether  point 
3  should  be  deleted  or  retained.  In  a  real  data  set 
many  points  could  be  expected  to  lie  near  the  critical 
value  of  point  3. 

The  line  shown  in  the  lower  plot  of  figure  6  cannot 
be  fitted  by  customary  least  squares.  It  can,  how- 
ever, be  routinely  fitted  by  the  iterative  procedures 
of  nonlinear  least  squares.  The  method  is  contained 
in  program  III  in  the  appendix.  The  coefficients  for 
natural  data,  under  the  restraint  that  predicted 
runoff  cannot  be  negative,  are  shown  in  the  third 
column  of  tables  6  and  7. 

Table  7  shows  that  for  Watkinsville  the  incorpo- 
ration of  the  runoff  restraint  tended  to  increase  the 
slopes  of  the  two  lines  (increase  in  coefficient  b  and 
in  coefficient  c)  and  to  decrease  significantly  the 
value  of  the  intercept  coefficient  a.  This  follows  the 
■results  expected  from  the  sketches  in  figure  6. 
Table  6  shows  that  this  single  restraint  did  not 
produce  satisfactory  results  for  Taylor  Creek.  Coef- 
ficient 6,  is  2.2952.  Physically,  this  value  may  not 
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INDEPENDENT    VARIATE  X 

Figure  6. — Effect  of  setting  negative  errors  to  zero. 

be  more  than  1.0;  otherwise  runoff  at  some  value  of 
rainfall  would  become  greater  than  the  rainfall.  The 
Taylor  Creek  data  were  rerun  under  the  additional 
restraint  that  the  coefficients  b  could  not  be  greater 
than  1.0.  The  results  in  table  6  under  this  fitting 
criterion  show  that  this  procedure  is  numerically 
feasible.  However,  the  change  in  the  coefficients  is 
quite  large.  The  distortion  is  due  to  some  extreme 
values  of  rainfall  occurring  in  the  short  record  of  11 
years.  It  is  expected  that  a  longer  record  would 
produce  better  averaging. 

The  last  column  in  table  6  shows  an  estimated  set 
of  coefficients  used  for  simulating  a  runoff  record 
for  Taylor  Creek.  These  will  be  treated  further  in 
a  later  section. 


SEASONALLY  CONTINUOUS 
STATISTICAL  FREQUENCY 
DISTRIBUTIONS 

In  order  to  use  equation  6  as  the  predictor  model 
for  generating  a  series  of  future  runoff  events,  it 
is  necessary  first  to  generate  values  of  synthetic 
monthly  rainfall.  The  usual  procedure  is  to  generate 
random  numbers  and  then  scale  and  transform 
these  numbers  to  represent  stochastic  monthly 


rainfalls.  Snyder  and  Wallace9  developed  the 
three-parameter  log-normal  distribution  as  a  func- 
tional variate  transformation  of  an  embedded  nor- 
mal distribution.  The  variance  of  the  embedded 
normal  distribution  was  standardized  at  unity.  The 
mean  of  the  embedded  normal  and  two  parameters 
in  the  transformation  function  were  evaluated  by 
nonlinear  least  squares.  When  so  defined  and 
evaluated,  the  three-parameter  log-normal  dis- 
tribution is  a  good  device  for  generating  stochastic 
rainfall  values.  One  first  generates  pseudorandom 
normal  numbers  of  a  computer.  These  have  zero 
mean  and  unit  variance.  They  are  transformed  to 
variates  of  the  embedded  normal  distribution  by 
simple  addition  of  the  parametric  mean  of  this  em- 
bedded distribution.  These  variates  are  then  trans- 
formed by  the  variate  transformation  function, 
after  which  they  can  be  considered  synthetic  rainfall 
values. 

The  log-normal  probability-density  function  used 
by  Snyder  and  Wallace  is 

p(v)=[V2Uk  (v-o)]-1 

exp  |  -y2     In  (v-o)  2  y  (1o) 

The  variate  transformation  function  is 

In  (v~o)=kx.  (11) 

In  equations  10  and  11,  x  is  the  variate  of  the 
embedded  normal  distribution  of  unit  variance  and 
mean  m.  The  variate  v,  to  be  regarded  here  as  a 
synthetic  rainfall  value,  is  dependent  on  x  by  equa- 
tion 11,  utilizing  the  mathematical  parameters  o 
and  k.  Snyder  and  Wallace  evaluated  the  three 
parameters  o,  k,  and  m  by  fitting  equation  10  to  a 
histogram  of  a  rainfall  record  using  nonlinear  least 
squares. 

The  usual  procedure  in  evaluation  of  statistical 
frequency  functions  such  as  equation  10  is  analogous 
to  the  procedure  for  monthly  discrete  regression 
coefficients  discussed  earlier.  An  advance  in  proce- 
dure discussed  below  is  the  development  of  season- 
ally continuous  frequency  distributions  analogous 
to  the  seasonally  continuous  regression  coefficients. 

In  the  preceding  section,  values  of  the  coefficients 
a,  b,  and  c  were  made  seasonally  continuous  by 
fitting  an  interpolating  function.  The  same  device 
was  used  to  make  the  statistical  frequency  distribu- 
tions seasonally  continuous.  Parameters  o  and  k  in 
equation  10  were  defined  by  interpolating  func- 
tions, using  the  same  arrangement  as  in  table  5.  The 


9Cited  in  footnote  4. 
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parameter  m  was  specified  as  having  the  same  value 
for  all  months.  This  says  that  the  embedded  normal 
is  the  same  for  all  months.  Now  in  simulation,  one 
generates  stochastic  values  of  the  embedded  normal 
as  outlined  above,  and  then  reflects  these  to  the 
proper  calendar  month  using  appropriately  interpo- 
lated values  of  o  and  k  in  equation  11. 

An  increase  in  efficiency  of  utilization  of  histori- 
cal data  results  from  this  seasonally  continuous 
stochastic  model.  If  equation  10  were  fitted  to  the 
rainfall  record  for  each  calendar  month,  36  pa- 
rameters would  be  evaluated.  By  seasonal  inter- 
polation for  o  and  k,  only  seven  parameters  are 
evaluated,  roughly  a  5:1  gain  in  efficiency.  The 
numerical  effect,  as  with  the  seasonally  continuous 
regression  coefficients,  is  to  produce  better  averag- 
ing and  a  smoothing  of  irregularities. 

The  procedure  may  be  viewed  as  fitting  12  sta- 
tistical frequency  distributions  simultaneously 
through  nonlinear  least  squares. 

Program  IV  in  the  appendix  is  a  subprogram 
written  to  fit  the  12  frequency  distributions  simul- 
taneously. The  monthly  rainfall  data  for  Taylor 
Creek  and  Watkinsville  were  organized  into  histo- 
grams by  months  as  shown  in  figure  7.  The  seven 
parameters  of  the  seasonally  continuous  version  of 
equation  10  were  then  evaluated.  The  specific  dis- 
tributions for  each  month  were  calculated  and 
superimposed  on  the  historical  histograms  in  figure 
7.  The  seasonal  left-right  shifting  of  the  frequency 
functions  is  due  to  seasonal  continuity  of  parameter 
o.  The  seasonal  change  in  peakedness  is  due  to  sea- 
sonal continuity  of  parameter  k. 

The  reader  should  fully  understand  the  fitting 
objective  in  figure  7.  It  will  be  noted  that  the  histo- 
grams have  a  base  of  16  classes.  In  fitting,  class 
errors  are  calculated  for  12x16=192  classes.  The 
sum  of  squares  is  minimized  for  the  192  classes. 
Figure  7  does  not  illustrate  a  best  fit  for  any  particu- 
lar month,  but  an  overall  best  fit  of  12  months.  The 
16-class  base  was  used  because  this  was,  by  coinci- 
dence, the  width  required  for  both  watersheds,  for 
September  in  Taylor  Creek  and  for  July  and 
November  at  Watkinsville.  The  concept  of  empty 
classes  beyond  the  range  of  historical  data,10  thus 
applies  to  all  months  except  these  three. 

Numerical  values  of  the  parameters  derived  from 
the  historical  data  are  listed  in  table  9.  Two  sets  of 
parameters  are  shown  for  each  drainage  area.  The 
origin  of  the  weighting  values  in  the  table  is  given 
below. 

I0Snyder,  cited  in  footnote  4. 


Grant11  investigated  the  theoretical  basis  of  fit- 
ting distribution  functions  by  nonlinear  least 
squares.  Two  findings  are  important  here.  The 
method  is  identical  to  maximum-likelihood  estima- 
tion for  classified  data  if  the  class  errors  are  nor- 
mally distributed.  A  bias  in  estimation  can  be  re- 
duced by  using  weighted  nonlinear  least  squares. 

The  weighting  of  errors  used  by  Grant  is 

eu.c=[hc-n(v)(.y[n(v)cf .  (12) 

In  equation  12,  the  usual  statistical  error  is  the 
observed  number  in  the  class  hc,  minus  the  average 
value  of  the  calculated  distribution  function  for  the 
class  n( v)c.  This  error  is  weighted  by  division  by 
n(v)c  exponentially  scaled  by  the  parameter  y. 
Grant  pointed  out  that  if  y  is  zero,  the  weighted 
errors  for  a  class,  ewc,  become  identically  the  un- 
weighted errors.  If  y  is  unity,  the  method  of  non- 
linear least  squares  becomes  identical  with  the 
method  of  minimum  chi-square.  Grant  found  that  a 
value  for  y  of  0.75  minimized  the  bias  in  estimation 
of  the  parameters  in  the  two-parameter  gamma 
distribution.  The  values  of  y  best  suited  to  other 
distributions  is  apparently  completely  unexplored. 

The  two  sets  of  parameters  in  table  9  are  the 
results  of  fitting  by  weighted  nonlinear  least 
squares,  with  y  set  to  0.50  and  0.75.  The  two  sets  of 
parameters  for  the  two  watersheds  in  table  9  are 
different,  but  the  only  relatively  large  difference  is 
in  o5  for  Watkinsville.  If  equation  12  is  solved  for  v, 

nGrant,  J.  L.  1973.  Statistical  frequency  analysis  by  optimiza- 
tion of  density  functions  to  histograms.  176  pp.  Doctoral  disser- 
tation, Georgia  Institute  of  Technology,  Atlanta. 

TABLE  9. — Parameters  of  seasonally  cyclic 
distribution  functions 


and 


parameter  0.50  0.7_5 


Taylor  Creek: 

o, 

-1.5545 

-1.6005 

°5 

.7175 

.7300 

Oq 

.3378 

.4530 

fc, 

.4953 

.5250 

k-a 

.7976 

.7663 

ka 

.8097 

.8042 

m 

1.8066 

1.8176 

Watkinsville: 

.5970 

.5296 

os 

.3909 

.1981 

-.5786 

-.6066 

fc, 

.7124 

.6903 

k-, 

.7695 

.7440 

k9 

.6228 

.6130 

ni 

1.7534 

1.8159 
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RAINFALL    IN  INCHES 

FiGl'RK  7. — Seasonally  continuous  rainfall  distributions. 
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Figure  8. — Runoff  simulation  comparisons  for  Watkinsville  watershed. 
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the  variate  transformation  function  takes  the  form 
v—exp  (kMx)+oM.  (13) 

The  two  different  values  for  o5  in  table  9  thus  take 
on  the  physical  significance  that  a  simulated  rainfall 
value  is  changed  by  about  0.2  inch.  This  does  not 
seem  of  major  significance  for  a  stochastic  variate  in 
an  observed  range  from  0  to  16  inches.  It  is,  how- 
ever, an  example  of  the  detailed  exploration  neces- 
sary in  further  development  of  weighted  nonlinear 
least  squares. 

The  parameter  values  obtained  with  y =0.75  were 
used  in  the  following  section. 

GENERATION  OF  SYNTHETIC 
RUNOFF  SERIES 

Programs  I  through  IV,  discussed  in  previous 
sections,  are  essentially  analytical.  That  is,  they  are 
intended  for  analysis  of  historical  data,  to  extract 
information  in  the  form  of  numerical  values  of  co- 
efficients and  parameters.  Program  V  is  not  a 
data-analysis  program;  it  brings  together  elements 
from  the  preceding  programs  to  synthesize  new 
information. 

The  basic  synthesizing  mechanism  in  program  V 
has  equation  6  as  its  core.  However,  researchers 
can  easily  modify  this  program  if  it  appears  desir- 
able, after  processing  a  particular  historical  record 
through  programs  I  through  IV,  to  use  another  of 
generating  equations  4  through  9. 

The  full  generating  mechanism  for  synthesis  re- 
quires adding  an  error  to  equation  6,  as  in — 

R  0  u  =aM + bslPsl  +cROu  _ ,  +e .  (14) 

This  error  is  implicit,  of  course,  in  the  very  defini- 
tion of  least  squares.  It  is  the  error  residual  to 
fitting  a  set  of  equations  to  a  historical  data  set. 
Such  errors  are  usually  assumed  to  be  normally 
distributed. 

In  using  equation  14  to  generate  a  synthetic  data 
series,  both  e  and  rainfall  Pu  must  be  considered 
stochastic  variates.  They  are  calculated  from 
computer-generated  pseudorandom  numbers,  by 
the  following  method. 

Normally  distributed  variates  may  be  standard- 
ized to  zero  mean  and  unit  standard  deviation. 
Standardization  is  performed  by  subtracting  the 
mean  from  each  item  in  a  sample  and  dividing  the 
difference  by  the  standard  deviation.  Since  any  two 
standardized  variates  must  be  equal,  one  can  write 
immediately 

r-mr  =  e-me  _  (15) 
sr  se 


In  equation  15,  r  is  the  pseudorandom  number, 
mr  is  the  mean  of  r,  sr  is  the  standard  deviation  ofr, 
e  is  the  stochastic  variate  of  equation  14,  me  is  the 
mean  of  e,  and  se  is  the  standard  deviation  of  e.  If  mr 
is  zero  and  sr  is  unity,  then  e  is  expressed 

e=rse+me.  (16) 

Equation  16  converts  a  computer-generated 
pseudorandom  number r  to  the  stochastic  variate e. 

Theoretically,  me  is  zero  when  the  method  of  least 
squares  is  employed.  For  nonlinear  least  squares  me 
should  also  be  zero.  In  practice  it  is  sometimes  very 
small,  but  not  zero.  This  nonzero  me  is  the  differ- 
ence between  the  means  of  calculated  values  and 
means  of  observed  values  in  tables  6  and  7. 

The  rainfall  stochastic  variate  is  generated  from 
the  continuous  seasonal  distribution  of  the  preced- 
ing section.  It  should  be  recalled  that  the  distribu- 
tion is  based  on  an  embedded  normal  distribution. 
The  standardizing  equation  is 

r—mr  _x—mx 

Now,  if  mr  is  zero  and  sr  is  unity  as  before,  andSj., 
the  standard  deviation  of  the  embedded  normal,  is 
unity  by  definition,  then 

x=r+mx.  (18) 

Substituting  equation  18  into  the  variate  trans- 
formation equation  and  specifically  using  P  for  pre- 
cipitation in  place  off  gives 

P,,  =  exp  [kJr+m)]+ou.  (19) 

With  e  and  Pv/  generated  from  pseudorandom 
numbers  r,  using  equations  16  and  19,  equation  14  is 
readily  evaluated. 

Figure  8  illustrates  the  results  of  runoff  simula- 
tion for  Watkinsville.  The  monthly  means  of  gener- 
ated data  are  shown  at  the  top,  and  standard  devia- 
tions are  shown  at  the  bottom.  To  the  left  are  results 
based  on  coefficients  derived  by  using  natural  data, 
unrestrained,  taken  from  table  7.  To  the  right  are 
results  based  on  coefficients  derived  by  using 
natural  data  with  calculated  negative  runoff  values 
set  to  zero  during  fitting,  also  from  table  7.  The 
positive  bias  in  synthetic  data,  illustrated  by  the 
upper  plot  in  figure  6,  is  clearly  shown  in  part  A  of 
figure  8.  The  generated  means  of  50  years  for  four 
sets  lie  above  the  means  of  record.  In  part  C  of 
figure  8  are  the  means  of  10  sets  of  50  years  each, 
using  coefficients  resulting  from  restrained  fitting. 
Here  the  means  have  been  lowered  towTard  the 
mean  of  record. 

It  is  not  possible  to  see  the  same  effect  in  the 
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standard  deviations  in  parts  B  and  D  of  figure  8.  The 
standard  deviations  were  lowered  for  the  months 
of  August,  September,  and  October,  but  not  appre- 
ciably for  other  months. 

Fiering  and  Jackson12  state  that  first  and  second 
moments  should  be  preserved  during  simulation. 
This  requirement  is  not  met  in  figure  8,  possibly 
because  of  the  relative  lengths  of  real  and  simulated 
records.  Fontane13  found  somewhat  similar  results 
with  respect  to  sampling  the  third  moment  of  a 
gamma  distribution. 

The  effect  of  record  length  on  distributions 
bounded  on  one  side  may  be  nonrigorously  de- 
scribed as  follows.  Consider  an  infinitely  long  rain- 
fall record.  Within  this  record,  the  items  are 
bounded  by  a  low  value  of  zero.  Scattered  through 
the  record  are  a  few  very  large  values.  Now  con- 
sider a  short  time  span  and  a  long  time  span  to  be 
selected  from  the  infinitely  long  record.  The  chance 


I2Cited  in  footnote  2. 

13Fontane.  D.  G.  1970.  Statistical  tolerance  limits  for  a  Pear- 
son type  III  distribution.  104  pp.  Master's  thesis,  Georgia  Insti- 
tute of  Technology.  Atlanta. 


of  including  a  very  large  event  is  greater  for  the  long 
sample  than  for  the  short  sample.  The  same  reason- 
ing indicates  that  a  long  simulated  record  is  more 
likely  to  contain  a  very  large  event  than  a  short  real 
record. 

If  the  distributions  were  unbounded,  one  could 
expect  extreme  positive  values  to  balance  extreme 
negative  values  as  sample  size  increased.  No  such 
balancing  is  possible  with  a  distribution  bounded  on 
one  side.  The  effect  is  shown  graphically  in  figure  9. 

In  figure  9,  histograms  of  rainfall  are  compared 
for  three  different  calendar  months.  The  observed 
record  length  is  25  years.  The  simulated  record 
length  is  50  years.  In  February,  the  observed  high- 
est monthly  rainfall  was  between  9  and  10  inches. 
For  the  simulated  record  the  highest  value  doubled, 
to  20-21  inches.  For  June  the  highest  simulated 
value  is  also  nearly  double  the  observed  value. 
These  extreme  values  exert  tremendous  influence 
on  the  moments  of  the  samples.  How  much  of  this 
effect  is  present  in  figure  8  cannot  be  answered  at 
this  time,  because  it  cannot  be  known  whether  the 
seasonal  log-normal  distribution  is  in  fact  the  proper 
distribution. 


TABLE  10.— Results  of  simulation  for  Taylor  Creek  watershed 


Month 


Observed 


200-year  simulation 


(2) 


(3) 


(4) 


Jan                                       0.451  0.786 

Feb  526  724 

Mar  926  .786 

APr  220  1.030 

May  360  1.428 

June                                      1.655  1.467 

My                                      1.679  2.676 

AuS                                      2.072  2.991 

SePl                                     3.210  2.936 

0ct                                       2.157  2.810 

Nov  308  1.584 

Dec   A62  1.031 

lan                                       0.725  0.906 

Feb  828  944 

Mar                                      1.350  1  124 

Apr  270  1.143 

^  650  1.760 

JrU™                                      2.719  1.820 

JAUly                                      1.807  3.065 

Aug                                    1-985  3.864 

XCpt 3-384  3.716 

°,Ct                                       2.971  3.235 

*0V  539  1.943 

Dec  135   1.193 

'Natural  data— unrestrained. 

2RO^Q  required  during  coefficient  evaluation. 

3«0&Oand  frssl  required  during  coefficient  evaluation 

Estimated  coefficients  (table  6). 


Mean  runoff 


0.506 
.268 
.083 
.106 
1.182 
1.521 
1.971 
2.614 
3.340 
2.840 
2.440 
1.192 


0.008 

0 
0 
0 

.130 
2.774 
4.086 
4.426 
3.771 
1.699 
.381 
.058 


0.101 
.055 
.311 
.895 
1.466 
1.842 
2.710 
3.006 
3.091 
2.710 
1.342 
.952 


Standard  deviation  of  runoff 


3.148 
2.146 
.986 
.687 
3.083 
2.670 
3.241 
3.192 
4.041 
4.207 
5.091 
4.125 


0.117 

0 
0 
0 

1.467 
4.253 
4.479 
4.064 
4.212 
3.441 
1.812 
.814 


0.484 
.255 
1.101 
1.729 
1.896 
2.245 
4.208 
3.575 
4.639 
3.949 
2.193 
.398 
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The  uncertainty  described  above  is  caused  in  part 
by  using  one  short  record.  Similar  uncertainty  re- 
sulted when  runoff  was  simulated  for  Taylor  Creek. 


1 0 1— 


0         2        4       6       8       10      12      14       16       18      20      22      24  26 


INCHES 

Figure  9. — Simulated  and  recorded  rainfall  distributions  for 
Watkinsville  watershed. 


The  resulting  average  values  of  means  and  standard 
deviations  are  shown  in  table  10.  Four  different 
200-year  simulations  were  performed  correspond- 
ing to  the  four  sets  of  coefficients  derived  from 
untransformed  data  and  listed  in  table  6. 

Not  one  of  the  four  simulation  results  in  table  10  is 
very  close  to  the  record  averages  of  means  and 
standard  deviations,  but  the  record  length  of  real 
data  was  only  11  years.  One  is  hard  pressed  to 
decide  whether  the  simulated  values  are  rational 
future  expectations.  Within  the  criteria  specified 
they  are  valid  statistical  expectations. 


CONCLUSIONS 

Synthetic  data  such  as  generated  by  this  package 
of  five  computer  programs  for  analyzing  and  syn- 
thesizing time  series  data  are  needed  for  use  with 
watershed  models,  to  predict  possible  runoff  of 
water,  sediment,  and  chemicals. 

The  illustrations  do  not  represent  a  fully  tested 
set  of  procedures,  but  are  intended,  rather,  to  show 
the  difficulties  met  when  a  synthesis  is  based  on  the 
very  short  record  available  for  the  average  research 
watershed. 

We  hope  that  other  researchers  will  be  interested 
in  applying  these  procedures,  or  modified  versions 
of  them,  to  simulation  problems  on  other  water- 
sheds. Pooling  of  results  should  make  possible  much 
better  tests  and  refinements  of  procedures  than 
could  be  obtained  with  the  limited  resources  at  one 
location. 


APPENDIX.— COMPUTER  PROGRAMS 


Program  I 
Input  variable  definitions 


NSET 
N 

LAG 


RO 


PR 


Number  of  data  sets  to  be  run. 
Number  of  observations  in  a  data 
set. 

Number  of  lags  for  which  auto-  and 
cross-correlation  coefficients  are 
to  be  computed.  LAG=1  is  zero 
lag. 

Values  of  the  dependent  variable 
(runoff). 

Values  of  the  independent  variable 
(precipitation). 


Output  variable  definitions 

SROl  Sum  of  the  lagged  dependent  values. 

SSR1  Sum  of  squares  of  SROl. 

SR02  Sum  of  the  unlagged  dependent 

values. 

SSR2  Sum  of  squares  of  SR02. 

SPR1  Sum  of  the  lagged  independent 

values. 

SSP1  Sum  of  squares  of  SPR1. 

SPR2  Sum  of  the  unlagged  independent 

values. 

SSP2  Sum  of  squares  of  SPR2. 

SSRO  Sum  of  products  of  SROl  with  SR02. 


19 


SSPR  Sum  of  products  of  SPRl^ithSPR2.  CORPR         Autocorrelation  coefficients  for  in- 

SPPR  Sum  of  products  of  SROl  with  SPR2.  dependent  variable. 

CORRO        Autocorrelation  coefficients  for  de-  CPRRO         Coefficients  of  cross  correlation  of 
pendent  variable.  variables. 


C  COMPUTATION   OF  AUTO-  AND  CROSS-CORRELATION  COEFFI^IEMTS 

OIMENSI  ON  R0(  1000)  ,  PR  1 1 000  )  ,  S  R01  (  40  ) ,  SR02  ( 43  ) ,  SSR  1  ( 40  J ,  S  SR2 (  40  J  ,  SP 

1RH40),  SPR2(40)  ,SSP1(40)  ,  SSP  2  (  40 )  ,  SS=*  Q  (  4  0)  ,  SSP  R  I  4  0 ) ,  S  PPR  (40  )  ,  CORRO 
1 (40) , COR PR (40) , CPRRO (40) 


00  111   LS  ET  =  1 , NSET 

READlr-,100)    N  *  LAG 

100 

FORMAT (215) 

READ  I  5  » 1 01 )    (R0(I)  ,1=1, N) 

KE  AO  ( i> »  101)  (PK(  I)  ,  1=1, N) 

101 

F0RMATI20F4.2) 

00   lG<t   1  =  1, LAG 

SROl ( I  )  =  Q.G 

SR02 ( I ) =o.O 

SSRK  I  )  =  0  .  0 

SSR2( I) =0.0 

SPR1(I)=0.0 

SPR2 (  I  ) =0 .0 

SSP 1 ( I )=0.0 

SSP2( I ) =0.0 

SSRGl I i =0.0 

SSPRI I)=0.0 

144 

SPPRt I )=0.0 

DO   106    1=1, LAG 

NEND=N- 1+1 

DO  106  J=1,NEND 

SROl I  I) =SR01 ( I ) +RG ( I +  J-1 ) 

SSRK  I  )  =  SSR1(  I  )+R0( I+J-l) *Ru(  I  +  J-l) 

SR02I I)=SR02( I ) +R0 ( J ) 

SSR2( I )=SSR2( I ) *RO( J)*R0( J) 

SPR1 1 1)=SPR1 ( I )+PR( I+J-l) 

SSPK  I)=SSP1(  I)+PR(I+J-1)*Pk(  I  +  J-l) 

SPR2 ( I)=SPR2( I ) +PR ( J ) 

SSP2 1 1 ) =SSP2( I ) +PR ( J ) *PR( J) 

SSRG( I )=SSR0( I )+RGl I+J-1)*R0( J) 


SSPR ( I ) =SSPR{ I)+PR (I+J-l ) *PR ( J ) 
10b  SPPR  (  I  )  =  SPPM  I )  +RG  ( I+J-l) *PR(  J) 

WRITE(6,112) 
112   FORMAT ( *Q   SUMS  AND   SUMS  OF   PRODUCTS* ) 
WRITE(6,108)    I I,SK01( I ) ,1=1, LAG) 
WRITt(o,l06)    ( I. SSRK I) .1=1. LAG) 
*RITt(b,l08)    I  I , SR02 ( I )  , I  =  1 , LAG) 
wklTE(6,108)    ( I,SSR2 ( I ) , 1=1, LAo) 
WRITE(6,108)    (  I.SPRKI)  ,1=1, LAG) 
WRITE (6, 108)    ( I,SSP1 ( 1 ) , 1=1, LAG) 
WRITE(6,108)    ( I ,SPR2 ( I ) ,1=1 , LAG) 
*KITE(6,108)    ( I ,SSP2( I ) ,1=1, LAG) 
WRITE(6,108)    ( I,SSR0(I ) ,1=1, LAG) 
WRITE(6,108)    ( I ,SSPR( I ) ,1=1, LAG) 
WRITE(6,108)    ( I,SPPR(I),I=1,LAG) 
108  FURMAT ( *0* ,   61 15.E16.8) ) 
DO   107   1=1, LAG 


20 


,  NMI=N-I+1 

CORRO ( I )  =  ( SSRO( I )- ( 1.0/ (NM1 ) ) *SR01 ( I ) *  SR02 ( I )  )  /  SwRT(  ( SSR2 ( I ) - ( 1 .  0/ 
1  (NMI  )  )*SRG2 ( I )*SR02(  I )  ) *{ S SR 1 ( I ) - ( 1. 0/ I  MM  I )  1*SR01 ( I ) * SR01 1 1 ) )  ) 

CORPRl I  ) =  ( SSPRt I )-( 1 .0/ (NMI ) ) *SPR1 ( I )*SPR2( I)  ) /SJKT(  ( SSP21 I  )- (  1. 0/ 
HNMI  )  )*SPR2l  I  )  *SPR2(  I ) )*( SSP1 ( I )-( 1. 0/ (NMI ) ) *SPR1 I  I) *SPR1 (I ) ) ) 

CPRRG{I)=( SPPR(  I)-{1  .0/  (NMI)  )  *SPR2(  I )  *SROi(  I  I  )  /S«iR  T(  (  SSP21I  )-(  1.0/ 
1  (NMI)  )  *SPR2 ( I ) *SPR2( I  )  )  *  (  S  SR  1  ( I  )  -  {  1 .  0/  { NMI  )  )  * S2  CI  I  I) *SR01 (I ) ) ) 
107  CONTINUE 

WKITE(6, 109) 

109  FORMAT ( *0     AJTO-CORRELATI ON  OF  DEPENDENT  VARIABLE*) 
*KITE (o, 105) ( I » CORRO ( I ) ,1=1, LAG) 

105   FQRMAT(*0*»10( I  5, F7.4 ) ) 
WRITE(t>  ,110) 

110  FGRMATl.*0  AUTO-CORRELATION  OF    INDEPENDENT  VARIABLE*) 
WRITE (6 , 105) (  I , COR  PR (I), 1=1, L AG) 

WRITE(t»,  113) 
113   FORMAT ( *0  CROSS-CORRELATION  OF  VARIABLES*) 
*RITE(6,105)( I,CPRRO(I) ,1=1, LAG) 

111  CONTINUE 
END 


Program  II 
Input  variable  definitions 


NSETS         Number  of  data  sets  to  be  run. 
NPTS  Number  of  observations  in  a  data 

set. 

RO  Values  of  the  dependent  variable 

(runoff). 

PR  Values  of  the  independent  variable. 

Output  variable  definitions 

CASE  A 

Al  SumoftfCVj. 

A2  Sum  of  ROM. 

A3  Sum  of  squares  of  POu _,. 

A4  Sum  of  products  of  POw_j  with  ROM. 

A5  Sum  of  squares  of  ROM. 

ERROR        Observed  minus  predicted  RO^. 

AA  Coefficient  a . 

BB  Coefficient  c. 

CASE  B 
Al  Sum  of  PR  u. 

A2  Sum  of  PO.u. 

A3  Sum  of  squares  of  PRM. 

A4  Sum  of  products  of  PR  M  with  ROM . 

A5  Sum  of  squares  of  ROM. 


ERROR  Observed  minus  predicted  RO^,. 

AA  Coefficient  a. 

BB  Coefficients. 

CASEC 

Al  Sum  of  Pi?  w. 

A2  Sum  of  PO',,^. 

A3  SumofPGV 

A4  Sum  of  squares  of  PP.V. 

A5  Sum  of  products  of  PR  M  with  POA/. 

A6  Sum  of  products  of  PPV/  with  PO>, 

A7  Sum  of  squares  of  POA/_  , . 

A8  Sum  of  products  of  ROM^  withPO 

A9  Sum  of  squares  of  P0V. 

ERROR  Observed  minus  predicted  PO(/. 

AA  Coefficient  a. 

BB  Coefficients. 

CC  Coefficient  c. 


Note  1. — On  second  pass,  values  for  cases  A,  B, 
and  C  are  based  on  log-transformation  of  PR  and 
RO. 

Note  2. — For  case  A  and  case  C,  first  output  set  is 
for  February,  and  last  output  set  is  for  January.  For 
case  B,  first  output  set  is  for  January,  and  last 
output  set  is  for  December. 


C   CALCULATION   OF  AUTO-  AND  CROSS-REGRESSION 
C  COEFFICIENTS  BY  MONTHS 

DIMENSI  ON  R0(  100 J)  ,_PR 1 1  000  )  ,  A 1 (  1 2  )  ,  A2  (  12  ) ,  A3  (  I  2  ) ,  AM  1  2 )  ,  AM..12)  ,  Ao( 
112) ,A7( 12) ,A3( 12) , A9( 12) .ERROR (100) 

READ(5,1000)  NSETS 

DO   1053   LSET=1, NSETS 

READ(5,1000)  NPTS 
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Iu04 
1003 

1007 


100t> 
1005 


1000  FORMAT (14) 
RtAD(5,1001)    ( RO ( I ) » 1  =  1 » NP  TSJ  

1001  FURMAT ( 20F4.2 ) 

READ  I 5t 1001  )    I PR( I ) , 1  =  1 ,NPTS)  

ASSIGN   1020   TO  KSPOT 

N=NPTS-1   

C   CASc   A.    bACK  RUNOFF  UNLY 

1022   DO  1002  11=1,12   

AH  1 1 1  =  6.6 

A2(II)=0.0   

A3( I  I  )=0.0 

A4(II)=0.0   .  

1002  A5(II)=0.0 

DO   1003   11=1,  12    5  

DO   10G<*    1=11, N, 12 

All  I  I)  =  Al ( I  I ) +R0{ I  ) 

A2  (  I  I )  =  A2 (  II) +R0{ I +1 ) 

A3(  I i )  =  A3 I  I  I )+R0( I )*R0{ I } 

A4( I  I )  =  A4( I  I) +R0<  H-1)*R0(  I ) 

AM  I  I  )  =  A5  (  I  Ii+RO{  1+11 *RQ(  1+1)  

CONT INUE 
WR I T  E ( 6 ,  1007 ) 

FORMAT ( *0  VALUES  FOR  CASE  A*) 
DO   1005  11=1,12 

WRITE {6, 1006 J    I  1, Al (  I  I } , A 2 ( I  I }, A3 ( I  I ) , A^ I  I  I ) , A5( I  I) 

FOR  MAT ( *0*, I3,5E17.9) 

CONT I NU  E 
vO   1 028  11=1,12 
IF( I  l.EQ.12)    GO  TO  1026 
NUM=NHTS/12 
GO   TO  1027 

1026  NUM=  (NPTS/12J-1   

1027  X1=A1 { I  I ) /NUM 
X2=A3 (I  I  )/ A  1  (  II  ) 
X3=A2 { I IJ/NUM 

X4=A4( I  I )/Al( I  I  )   

bb= ( X3-X4)/ (X1-X2 ) 

AA=X3-l BB*X1)     

WRITE (6  ,  1050) 
1050  FORMAT ( *  0* , / ) 
IX=0 

DO   1030  IERR=II,N,12 
I  X=  I  X+  1 

103U   ERROR!  IX)     =RQ(  IERR+1)-AA-(BB*R0<  i ERR  )  )  

WRITEi6, 10311    I  I , (  I  ERR, ERROR  I  I EKR ) , I  F  R  R= 1 » N  J  M ) 

1031   FORMAT! *OERRORS  FOR* , I  5/ ( 8 { 1 4  ,F 10 . 7)JJ_ 
WRI TE(6, 1029)  AA,BB 

1029   FORM AT ( *0  AA    IS*,F16.8,*       Bd  IS*,F16.6) 

lU2o  CONTINUE 
C   CASE   B.   CURRENT   PRECIPITATION  ONLY 
DO   1008  11=1,12 

A1(U)  =  0.0   

A2( I  I ) =0. 0 
A3  (  1 1 )  =  0  •  0 
A4(  I  I )  =  0.0 

lu08  A5(II)  =  0.0  

DO   1009    11=1,  12 

DO   1010  I=II,NPTS,12 
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AH  I  I  )  =  A1(  I  I)  +PR(  I  ) 


A2( I 1)=A2( II)+RO(  I  J 

A3( I  I  J  =  A 3 { I  I )  +PR( I  l*PR< I ) 

A4  (  I  I  )  =  A4  (  I  I  )  +  PR(  I  )  *RGl  I  ) 


10  10 

I  U  i  w 

Ad(  I  I  )  =  A5  (  I  I  ) +R0(  I  )  *R0  (  I  J 

1U09 

GUN T  I NL)  F 

WR I TE (6  . 101 1 ) 

n  r\  4  •  L  \  u  y  a  w  a  a  / 

1  O  1  1 

F(ikN<AT(*0    VAI  UF^   F  f)R    fA^F    fi*  1 

DO   1012  11=1,12 

WKlTE(6tl006)  II,A1(II}.A2(1I)»A3(II),A4(II),A5(II) 

»  »  i  »  A    i    i—.*  vy  y  x  vy  \y  w  §         a    l    y        a.   v    j.    a   /    y  /~\  < —    \    i.    i    /   y  r-%        i    i.   a    /    y   '  *   ■    \    i    i.    /    y   ^  _y  %   a     i  i 

1012 

CONT  INU F 

DO   1 032    1  I = 1 •  1  2 

L •  —        J-   o  —  i—         XI        X.   f    X.  C. 

NUM=NPT  S/ 12 

A  1.      MX  \  X   1  J  /  INU1  1 

X2  =  A3 ( I  I  )  /Al (  I  I  ) 

X 3=  A2 (  I  I  )  /NUM 

X4=A<*  (  I  I  }  /Al  (  1  I  ) 

/A    l         HT  \    1    X    /   /    r-\  X    h    X    A  y 

BB=(X3-X4)/(X1— X2) 

AA=X3-{ 8B*X1 ) 

wR ITFIh.1050) 

IX=0 

DO   1033    I  ERR=  I  I  .  NPTS ■ 1 2 

I  X=  I  X+  1 

1033 

ERROR ( I X  )     =  R0(  1ERR)-AA-{B8*PP(  IERR)) 

wRl T  E ( 6  , 1 03 1 J    I  I , (  1  ERR , ERROR ( I ERk ) ,  I ERR= 1 , NJM J 

•*RITF(ft.  10791    A  A  .  ft  R 

1 03  2 

CONT  INUE 

C  CASE 

C.    CURRENT   PRECIPITATION   AND   BACK  RUNOFF 

UO  1013    1 1=1. 12 

A  1  (  I  I  )  =  0  .  0 

^  XI    AAV        <J  »  \J 

A2  (  1  I  )  =  0  .0 

A3! I T 1=0-0 

A4(  I  I )  =  0. 0 

A5 (  I  11  =  0.0 

Ac(  I  I  1  =  0.0 

A7I I I 1=0.0 

"l%XX7        \y  •  Vy 

A6 {  I  I  )  =  0 -0 

*-y    *    X    A    m         w  •  \S 

1013 

A9(  I  I  )  =  0. 0 

y  i    a    a    w        vy  •  w 

00    1014    I  1=1,  12 

U0    1013  I=II,N,12 

Al (  I  1 )  =  A1(  1  I ) +PR(  I +1  ) 

A2( I I)=A2 (II) +R0( I ) 

A3(  I  I )  =  A3{ I  I ) +R0( I  +  1 ) 

A4( I I)  =  A4( I  I ) +PR( I  +  1J*PR(  I+iJ  

A51 I  I J  =  A5( I  I ) +PR( I +  1  )*R0(  I  ) 
A6U  I  )=A6(  II  )+PR(  I +1)  *RC(  1  +  1) 
A7( I  I  J  =  A7 ( I  I ) +R0( I )*R0( I ) 
■   Ao(  I  I )  =  A8( I I)+R0( I )*R0{ 1  +  1  ) 

1015  A9( Ili  =  A9(I I)+Rd{ I*1J*R01 1+1 J 
1014  CONTINUE 

WRITE(6,1016) 

1016  F0KHAT(*0  VALUES  F0k  CASE   C* ) 
DO   1025    1 1=  1,  12 

WRI  TE(6  ,  1017)    I  I  ,  Al  {  I  I  )  ,  A2  (  1  I  ),  A3  {  I  I)  ,  A4  (  I  I  )  ,  AM  I  I  ),  A6(  I  I  )  ,  A7(  I  I  )  , 
1A81I I), A9(II) 

1017  F0RMAT(*0*,  13  ,  ( 5E17.^  )  ) 
1025  CONTINUE 

U0   1034    11=1,  12 
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IF( I I.6Q.12)   GO  TO  1035 

NUM=NPTS/12 

GO  TO  1036 
1035  NUM=i NPTS/12J-1 
luJo  X1  =  A1 ( I  I )/NUM 

X2  =  A2UI)/NUM 

X3=A3( I  I )/NUM 

X4=A4( I  D/Ali  Ii  ) 

X5=A5 {  I  I )  /AH  I  I  J 

X6  =  At> ( I  I ) / Al (  I  I  ) 

X7=A5(  I  IJ/A2(  II  ) 

XS=A7( I  I  )/A2(  I  I ) 

X9=A8I  I  I  )/A2(  Hi 

Yl={ X2-X5)/{ X1-X4) 


Y2=( X3-X6)/ (X1-X4) 

V3=( X5-X8J/( X4-X7) 

■ 

Y4=(X6-X9)/{X4-X7) 

CC=l Y2-Y4)/(Y1-Y3 ) 

BB=Y2-iCC*Yl) 

AA=X3-(CC*X2 )-{ BB*X1 )  

WRI TE ( 6  t 1050) 

IX=0   

DO  103d  IERR=II,N,12 

IX=IX+1 

1U38 

ERROR (IX)  =R0(IERR+1)-AA-(6B*PR(IER^+1))-(C:*R0(IER<)) 

wR  I  T  E  (  6  ,  1 03 1 )    1 1 , ( IERR, ERROR ( IERR ) , I SRR= 1, NJM)  

WR I T  E { 6  t  1037)  AA,BB,CC 

103  7 

FORMAT ( *0  AA   IS*, Fib. 8,*       BB    IS*. Fib. 3,*       CC    IS*, Flo. 8) 

1034 

CONT  INUE 

AUGMENT   TO  PREVENT   LOG  OF   ZERO  ERROR 

GO  TO  KSPOT.I 1020, 1021) 

1020 

ASSIGN   1021   TO  KSPOT   

00  1018  I=1,MPTS 

PR  I  I  )= AL0G( PR (  I  J+0_*JL0_1 

1018 

ROU  )=ALCG ( RO I  D+0.01) 
HRITE(6, 1019) 

1019 

FORMAT ( *0  PRECIP  AND  RUNOFF    I  NCR E AS I J , 0 .  1 0  A^J    0.01,  Rc 

SPECTI VELY» 

1T0  PREVENT   LOG  OF    ZERO  ERROR*/*   FOLLOWING  ARE   BAStO  ON 

LOGS   DF  DAT 

2A*) 

WRI TE(6  ,1051 )    ( I KK, PR( IKK) , IKK=1 , NPTS ) 

1051 

FORM  AT ( *0  LOGS  OF   RAINFALL   PLUS   0. 10* / ( 8 { I  5 , Fl 0.6 ) ) ) 
rtRITE(6,1052i    (  IKK , RO ( IKK)  ,  IKK=l, NPTS) 

1052 

FOkMAT(*0  LOSS  OF   RUNOFF   PLUS   0. 0 1*/ ( 8 ( I  5 , F 1 0  .  o  )  )  ) 

GO  TO  1022 

1021 

CONT INUE 

1053 

CONT  INUE 

STOP 

END 
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Program  III 

This  is  a  subprogram  only,  intended  for  use  with 
some  other  program,  or  system  of  programs,  to 
obtain  optimum  values  of  the  coefficients.  The  pro- 
gram system  used  for  this  study  was  presented  in 
DeCoursey  and  Snyder.14 

To  remove  the  restraint  of  negative  values  set  to 
zero  during  fitting,  the  two  statements  labeled  A 
can  be  deleted. 

To  incorporate  the  restraint  of  coefficient  b  not 
larger  than  1.0,  the  short  LIMIT  subprogram  can 
be  used. 

To  fit  log-transformed  data,  the  four  statements 
immediately  following  the  program  END  can  be 
inserted  at  point  B. 
Input  variable  definitions 
OBS  Y  Values  of  runoff  by  months. 

RAIN  Values  of  rainfall  by  months. 


BRO  Back  runoff.  Runoff  in  the  month 

immediately  preceding  the  first 
OBSY  of  record  utilized.  It  can  be 
estimated. 

Other  input  controls  for  the  optimizing  procedure 
used  were  read  in  by  the  system  MAIN  program. 
Output  variable  definitions 

Output  is  controlled  by  the  system  MAIN  pro- 
gram or  other  subprograms.  Normally,  values  of 
the  coefficients,  calculated  values  and  errors  for 
each  observation,  and  sensitivity  coefficients  for 
each  observation  are  output  for  selected  rounds  of 
iteration  until  solution  converges  to  preset  iteration 
limit. 


14DeCoursey,  D.  G.,  and  Snyder,  W.  M.  1969.  Computer- 
oriented  method  of  optimizing  model  parameters.  J.  Hydrol. 
9(1):  34-56. 


ft  .  THlb  SUBPROGRAM  IS  USEJ  WITH   iTtKATIVE  NONLINEAR    LEAST  SQUARES 
C^TC   EVALUATE   SEAS  J N ALLY  CONTINUOUS  REGRESSION  COEFFICI ENTS. 
ft~         SUBROUTINE   PART  AL  t  NOBS »  N  PA P , A H U \ , K Pk 1 N , 1  ST o P ) 

(  COMMON/ ONE /PAR (  11 )  ,X(900,  12) ,CQEl II) ,0BSY(900)  ,PR EDY I  900 ) /TWO/C 3k  i 
1 12  1 12 )  , EVALll  1 )  jEVEClll.ll  ) 
K         DIMENSION  A13,12),   KAIN15 00)  ,RQl 50 0) ,PR(  1 1  ) 
jt~~.      IF(MRUN.NE.O)    GO  TO  1000 
M  HEAD (5,1001)    <  OBSY  1 1),  1=1  ,  MDBS) 

READ  lb, ^100U*RAIN{  I)  ,1=1, NOBS) 
_1L0L   FUKNipT  (20F4.2) 

READ<  5»  10C1 )   BRO.     B 

lGOu   Cu   ii.Ll    1=1  ,NPAR  

T"~       DR(I)=0.C  L*PAR  (  1) 

_  IF  I  A  H  S ( CR  t 1 ) )  .LT.O  .001  )    DR  I  I  1=0. JUl 

1CC2   cGi\T  I  N U E 
;jj  ASSIGN   2  008  TU    I  TR  AN  

rCGl   ASSIGN   1004  TO  KSPCT 

y  pi=pakmlj  

W~  >2  =  PAR(  2  ) 
Wt-  P3  =  PAM3) 
\~        1=  1 

}«007   Al  I  ,?.  )=  (-J.0937*P3  )  +  10. 8o  72*P  1)  +  (  Q.2266*P2) 
Wf  »     A I  I  ,  3)  =  1-0.  12  5*  Pi  )  *  J.  5  c  23*  IP1+P2 ) 

rW-        AU  ,4  )=  t-0  .0937*Pi:  )  *  i  J  .  I .1 t  u  v?  1  j  +  {J  .  t  o  /  2*  i  _  )   

AT1,5)  =  P2 

Al 1 , t)=(-0.0937*Pl )+(0. 6672*p^) + I C.22b6*P3l 
A(  i  ",7)  =  l-0.125*Pl  )+0\5625*TP2+pT) 
At  i  ,8)  =  (-Q.0937*Pl )  +  i0.226t*P2H-(l.cb72»pj) 
t  A(It9)=P3 
iMh       Al  I  ,l<J)  =  l-0.0937*P2)»(0.a6  72*P3)  +  10.2266*P1) 
^^Mpr^TTTTT)  =  (-0  .  1 2 5*P 2  )  +  0 .  5o~2 5*1  P  2+P  1 )  ~ 
\         All  fi2)-=  1-0.0937*  P2J_+J0  •  2  2 66_?  P3  )  *■  (0_.36 72  *P  I  J 
M  I  ,  1  )  =  P  1 

 GO  T  C   KSPCT  ,  (  1004,  10 05",  10  0  -  J 

1C04    ASiloK    UJ-)   TC  KSPGT 
P1=PAR(4) 
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P2=P Ak ( 5 ) 


SUBKLUT  INE  L  I  M  I T  ( M VAR  t MfcUN  )  ,   .   ■ "  % 

COMMON/ ONE  /PAR  (  1 1  A?  9  X  (   900 , 12 )  ,COE ( 1 1) , GBSY {  900),PREDYl   9  00) 

_  C°  1000    1=4 1 6  __ 

1F(PARI I ).GT.l .0)   GO  10  1001 

 GO  TO   1 0 U 0  .  

1C01  PAR  i  lT-1  .0  •  .,> 

 C0£(  I )  =  1 « O-PAR  1 1  )   _ 

1UU0   CONT  1NUF 


Program  IV 

This  is  a  subprogram  only,  intended  for  use  with 
the  same  systems  of  programs  as  program  III. 
Input  variable  definitions 

NYR  Number  of  years  in  the  record. 

NCLAS         Number  of  classes  for  each  month. 
WT  Value  of  exponential  parameter  for 


weighting  errors. 

OBSY  The  number  of  observations  in  each 

class.  Observations  are  read  in  as  a 
continuous  set  and  then  organized 
into  months  by  Do-Loop  3001. 

Output  variable  definitions 
Same  as  for  program  III. 


C  THIS   SUBPROGRAM  IS   UiED   WITH   ITERATIVE   NONLINEAR   LEAST  SQUARES 
G  TO  EVALUATE   SEASONAL.Y   CONTINUOUS   DISTRIBUTION  FUNCTIONS 
SUBROUTINE   PARTAi_(  NO  si  ,N  P  AR  ,  MR  UN  ,  <PRIN  ,ISTQP) 

C0iM0N/0NE/PAR(il),X(ojj,12),C0E(ii),0BSY(5jv),H'r;E0Y(3^j)/lrtJ/:.;( 
112,12)  ,EVAL(ll),E\/EC(ll,li) 

DIMENSION   DR (7) ,536(25,12) ,A<2,1?) ,FUN(1U) ,AVFUN(3G0) 

IF(  MRU  N«NE«0_)    GO   TO   3.  J 

READ (5 ,6000)  NYR,MCLAS»WT 
ojuj   FORMAT (2Ih,F6.  C) 

M03S=1 2*NCLAS 

DO  5101  J  =  l,12 

IE=J*NCLAS 

IB=IE-NCLAS+1 
5101   *EA 0(5, 310  0) (OBSY { I) ,1=19, IE) 
5iOC   FORMAT (25F3.Q) 

3IF=SQRT( 2*3. 1h15327) 

DO   3jv 1  1=1,12 

I B= (1-1 ) *NClAS+1 

DO  3JG2  IS=1,NClA5 
3jj2  SOddS  ,1)  =  03SY  (I3HS-1) 
3001  CONTINUE 
3u,i   ASSIGN   null    TO  ITRAN 

ASSIGN   3j18  TO  JTRAN 

DO  i+Obj  I=1,NPAR 

ORCI)=j.Q05*ABS(PAR(I ) ) 
'  IF(DR(I)  .LT.0.005)  DR<I)=Q.G05 
i»G60  CONTINUE 

DO   3ulJ   1=1, SjO 

0  B3  Y  (  I  )  = j  •  0 
3u lb  3  RE  DY ( I) =0  .  j 
3011  ASSIGN   30  0  3  TO  KS^OT 

P1=PAR(1) 

P2=PAR<2) 

P3=PftR(3) 

1=1 

5uj5  A  (1,1)  =Pi 

A (I ,2)=(-0.0937+P3)+(0.8672*Pl)+(0.2266*P2) 

A  (1,3) =(-0.l25*33)  +l.5625*(P1+°2) 

MI, h)^(-G.u9 37 *33)+(u. 22 66*Pl>+(G.8o7  2*P2) 

A  (I  ,5)  =P2 

A(I,b)=(-G.C937*3i)-t-(w.3672*D2)  +  U.22b6*P5> 
A(I,7)=(-u.l25*3l)+-. 562  5* (P2+P3) 
A(I,b)=(-u.u9.}7*3l)+(..22o6*P2)  +  (L.  667  2*P3) 

\  (1,9) -PS 

A(I,lG)=(-C.w9  3?*32)  +  (Ci67  2*P3)  +  (c.  22  6o*pl) 
4(I,11)=(-w.1l5*P2)+v.?d2  5*(P3+P1) 
& (I ,12) =(-0.G  937*:i2)+(G.  226  6*P3)  +  (L.6o72*?l) 
GO  TO   <SPOT,  (3  0T3,  3wC'4) 
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3 ...  3 

ASSIGN   3u  04  TO  <S30T 

31=PAR(4> 
P2=pAR<5) 

=»3=PAR(6) 
1=2 

3  ii  J  4 

GO  TO  iWS 
MGTOE=Q 

HO   3u  0  6  r10N  =  l ,  12 
I  ADV=0 

CL=(-0 .005) 
ICLAiS=J 

IPC A( 1 , MON) «LT  »CD    GO   TO  3107   

3d  u8 

GO  TO   5  u  j  1 

IF( A<1 ,MON) .LT.Cl)    GO   TO  3LC7 

3  j  g  7 

GO  TO  3jj9 
IAJ V  =  I AOV+1 

3  L=CL-  1.  Q 

GO  TO  3jG8   

P.  jl 

IF( A (1 ,MON) .GT . (3u+l.  u) )    GO  TO  5G0 

GO  TO  3009 

P  1]  u  - 

ICLASS=ICLASS+1 

IF(S03 (ICLASS ,M0M) .GT.fl)    GO  TO  3Cu 

Q 

I ADV=I A  DV -1 

CL=CL+1.j  

3  u  0  9 

GO  TO   5  0  3  1 
*ITOT=N  3_  AS  +  IAOV 

"UNO  =  0  .  J 

00  3012  I=l,NTOT 

JO   3ul3  l4TH=l,lj 
3l=Cl  +  J  .  It 

IF( CL • L£. A ( 1 >  M ON) )    GO   TO  3C15 

GO  TO  3ul6   

5  a  1  3 

FUN (I4TH) =0.0 

3j  lo 

FUN(I-fTH)=£XP(-J.p*((AL.33(CL-A(l,M 
FUN(ItTri) =  FUN (I%TH)/(PIF*A<2»MON)* 

ON) ) /  A  ( 2  ,  MON)  ) -PARC  7) >**2) 
(Cl-A  (1  ,MON))  )  

3.,  13 

C  O.MTI  MUt 
AVFUN(I)    =  O.G 

SIjI 

DO  6lQ 1  11=1,9 

AVPUN(I)    =   AVFUN(I)    +  FUN (II) /lQ.fl 

A\FUN(I)    =   AVFUm(I)    +    (FUNfc    +  FUN( 
A  VF  UN ( I )  =  A VFUN ( I )  *  NY  R 

lu))/2G,C 

<rUNC  =  F  JN  (  10  ) 

GO  TO   ITRAN,  (tt  Oil,  4u  24) 

*tj  11 

NGT  SAV  =  NG TOE 

DO   301^  <=l,NTOT 

3j  m 

IX=NGTOE+K 

PRE ( IX) =AVFUN(<) 

IF(  IAOV.vit.O  )    GO   TO  nig] 

GO  TO   +  101   . 

mo& 

I09=NGT0E  +  1  +  IA0\/ 
IO£=NGTOE+NTOT 

rc=i 

00  41G  2   10=106, IOE   ;    h.i  '  .  ::i 

DB5Y  (I  J)  =S08(IC  fWom 
IC=IG*1 

4101 

GO  TO  -1j3 
I 03=NGT0^+i 

I 0E=NGTOE+NT0T 
I C= IA  BS (I AOV) +1 

30   HlG-f   I0  =  I0BfI3E 
0BSYCI01 =S0Bf ICtMON) 

W  J  t 
^  1 G  5 

i  U-  1  L/+  1 

CONTINUE 

DO  <fiil  3  L  =  l  ,NTOT 
IX=NGTO£+L 

IF (1ST 3P. EQ.u)    G  0 
GO  TO  1031 

TO  1J0J 

1  o  J  U 

X (IX,8)=0BSY (IX)- 

3 REOY  (IX) 

1  u  0  1 

J    \J           \     <J          T  J   I  U 

X (IX, 3) =( OBSY  (IX) 

-PR£DY(  IX))/(AV/FUN(l)**WT) 

■tj  lo 

fi  D  M  T  T  M 1 J  - 

GO  TO  3006 

t     7  - 

DO   3. 1 7  M=T»NT0T 
I  X=  NG  T  Jc+  M 

)  j  1  / 

J  w  J  O 

% (IX, JCOtT={ avfun 
NGTOE=NGTOE+NTOT 

<  tt)  -PREDY ( IX) )/DR(JCOE) 

30  TO    JTRAN,  (.3018,3019) 

ASSIGN   -fu^1*  TO  ITRAN 
403S  =  fJoTO^ 

JCDE=i 

^A*  (JCOt)  =  PhR  ( J-iJ 
30  TO  3011 

I  )  +  D  ;  ( J3^>£) 

3  J  19 

PA*  ( JCOE) =SC0£ 
IF(  JC0£.  ;_Q.NPAR) 

SO   TO  30^1 

JC0£=J30- +1 
GO  TO 

5  j  •+ 1 

RETURN 
E  NO 

Program  V 

Program  V  MAIN  calls  subroutine  RANTAB. 
Subroutine  RANTAB  calls  subroutine  RNNG. 
Subroutine  RNNG  computes  one  random  normal 
number.  Subroutine  RANTAB  draws  at  random 
from  a  table  of  100  random  normal  numbers  with 
random  replacement.15  This  method  extends  to 
almost  inconceivable  length  the  finite  series  of 
pseudorandom  numbers  before  recycling  begins. 
The  integer  numbers  appearing  in  these  sub- 
routines are  for  a  CDC  6400  computer. 
Input  variable  definitions 


STNS 

TITLE 
NSETS 


The  number  of  station  simulations  to 
be  performed. 

Aphameric  problem  title. 

Number  of  sets  (synthetic  series) 
for  each  station  simulation  (for 
each  set  of  coefficients  and  param- 
eters). 


NYRS 

RN,N1,N2 
BRO 


CLASW 


SDRES 
BARERR 
COE 
PAR 


Number  of  years  in  each  simulated 
series. 

Beginning  random  numbers. 

Back  runoff.  Runoff  in  the  month 

immediately  preceding  first  month 

simulated. 
Class  widths  for  the  construction  of 

histograms  of  monthly  rainfall  and 

runoff. 

Standard  deviation  of  residuals,  Se. 
Average  of  residuals,  e. 

aV  &5>  &9>  ^5>  ^9>  ^i,  C5,  C9. 

Oj,  O2,  O3,  /c  j,  /c 2,  ^3,  yyi • 


15Grant,  cited  in  footnote  11. 


Output  variable  definitions 

SUMP  Average  of  simulated  monthly  rain- 

fall values. 

SUMP2         Standard  deviation  of  simulated 

monthly  rainfall. 
P  Values  of  simulated  monthly  rainfall . 

SUMR  Average  of  simulated  monthly  runoff 

values. 
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SUMR2         Standard  deviation  of  simulated  histograms  of  rainfall  and  runoff. 

monthly  runoff  values.  RN,N1,N2     Last  values  of  random  numbers  gen- 
RO               Values  of  simulated  monthly  runoff.  erated.  Can  be  input  for  subse- 

NINCL         Number  of  events  in  each  class  of  the  quent  runs. 


C   THIS   PRCCRAC   SIMULATED   MONTHLY   ROiNCFF    IN    INCHES   USING   SEASONALLY  CYCLIC 
C   COEFFICIENTS  OdTAINO   BY   METHOD  OF   CLNTI \UULS   PAR ABGLI C   I  N  TER  POL AT  ft)N« 
CI  MENS  I  CN  X(12),A(12),B(12),C(12),D(12),  E  ( 1 2  )  »  COE  i  9)  t  PA'R  i  7)  »N  I NCL  ( 
1100, 12) * P( 12, ICC) ,R0( 12 t ICQ)  ,  TI  TLE (  G)  , SUMP( l2)  ,  SU  VP2  (  12  f , SUMR { 12)  , 

2SCMR2(12)    

INTEGER    SET  ,  YR  ,  RN,STNS,ST!\ 

 READ(5,1042)  STNS  

1C42   FOkMAT ( I  A )  .  - 

STN=1    _   

1044   READ (5, 1G37  i (TITLE 1 1),  1  =  1,8) 

1037   FORM A T ( 8A 1 0 )   

WRITE (6,  10  38 ) ( T  ITLE I  I  )  , I-  1,8) 
1C38   FukMATI*!* ,30X,8A1G) 

READ ( 5,  10 CO)   N  S  E T S  , N  YR  S  , R N  , Nl  ,  N2 

1CCC  FORMAT  12 15  1 31 15.1   

REAG(5,1CC1)    BRO ,C  LASW , SORES ,  BARERR 
1001   FORMAT ( 10F6  .0  ) 

READ(5,1GC1)    (C0E(  I)  ,1=1  ,9) 

 REAC(5t  ICO'l)    (  P  AR  (  I )  ,1  =  1,7)   

fcRITE  {6  ,1039)    SDRES  ,  BRO  ,  RN  ,N  1  ,N2 
1C3S   FORM A  T ( *0  STD.   DEV.   OF   RANCCf   ELEMENT* »F 10. 6/*     BACK   RUNOFF  IS*,F5 
1.2,*    INChES.*/*     BEGINNING  P  AND CM   NOS.   ARE*, 3120) 
V>R  I  TE  (6  ,1040)  (  I  ,COE(  I)  ,1  =  1,9) 
104C"F0RMAT l*C  THE  COEFFICIENTS  ARE... */ ( 91  1 3  ,F 10 . o ) ) ) 

 feRITE(6, 1041) M,PAR {  I  )  ,  1=1,7)  

1C41   FORMA T( * C  THE   DISTRIBUTION   PARAMETERS   FOR   RAINFALL  ARE.  ■  •  */ (7  (  13, F 

llu.6)))   

C   CALCULATE  MONTHLY    VALUES   uF  COEFFICIENT  A. 

  AA=COE i  1) 

6B  =  CCEl 2  ) 

CC=CCEC3)   

ASSIGN    1C02   TO  KSHIFT 
1C05  X(1)=AA 


M2) 

=  (-L. 

C937*CC)+ (0.8c72*AA) +(0. 

22t>6*BB ) 

X  (3) 

=  (-0. 

12  5*CC  )  +  0.5625*( AA+BB) 

X  i  4 ) 

=  1-0. 

0937*CCI  +  ( 0 .2  266  * A A  j +(0. 

8672*BB ) 

X(  5  ) 

=  BB 

A  (6  ) 

=  1-0. 

0937 *AA)+ 10.86 72 *BB  )  +  lO. 

Z26t>*CC  ) 

XI  7) 

=  1-0. 

12  5*AA)+0. 5  62  5* ( 8B+CC) 

X  l&) 

=  (-0. 

0937*AA)  +  (0.2266*BB)  +  (0. 

8672*CC  ) 

X  9) 

=CC 

A(10)  =  l-C.093  7*BB)  *T0.8672*CC )  +  ( 0 . 22o6*A  A) 
X( 11 )=(-G.125*BB)+0.562  5* (CC+AA) 
X(  12)=(-0.093  7*BB)  +  ( G.  2266 *CC  )+ (  0.36  7  2*AA) 
GO  TC  KSHIFT,  (1002,100^,1CC7,  1009,1011) 

1002   CO   1003  1=1,12 

1C03   A( 1 )=X{ I ) 

ASSIGN   1004  TO  KSHIFT 
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C  CALCULATE  MONTHLY    VALUES   OF   COEFFICIENT   B . 

/iA=C0E(4) 

BB=CCE( 5  ) 

CC=C0E  (6) 

GO  TO  1005 
100 n  D 0  10  jo  I=i, 12 
IdO  b  3  (I  )  =X  (I) 

ASSIGN   1„Q7   TO  KSHIFT 
C  CALCULATE    MONTHLY   VA.JES   OF  COEFFICIENT  C. 

AA=COE (7) 

3B=GOE  (8) 

CC=COE (9) 

GO  TO  1005 
1j07  30   luu^  1=1,12 
1,j8  Z  (I)  =X  (I) 
C  CALCULATE   MONTHLY    VALUES   OF  PARAMETER  O. 

AA=PAR (1) 

3B=PAR  (2) 

CC=PAR(3) 

ASSIGN   1j09  TO  KSHIFT 

GO  TO  1005 
1009  OO   1010  1=1,12 
Julfl   0  (I )  =X  ( I) 

ASSIGN   1J11  TO  KSHIFT 
C  CALCULATc    10NTHLY    VA.JES   OF  PARAMETER  K. 

A  A=PAR  U) 

3B=PAR(5) 

:C=PAR(6) 

GO  TO  1005 
1-11   00   1012  1=1,12 

1012  E  (I )  =  X  ( I) 
SET  =  1 

C  SET  UP  RANDOM   NUMdER  TABLE 

C  THEN  CALCULATE  PRECI3.    AND  RUNOFF    BY  MONTHS. 
GAlL   RANTAB(RN,Ni , N2 , DRAW , 1) 
103A  DO   1013  YR=1, NYR5 
MON  =  l 

101**   3  AL  L  RANTAB(RN,Nl , N2 , DR4W , 2) 

P(M0N,YR)  =  (EXP(£(M0N)MDRAW  +  PAR(7)))+D(M0N))*ClASW 
IFJP(MON, YR) .LT.O. 0)    PCHON, YR) =0 . 0 
CALL  RANT  ABC RN , Nl , N2 »  DRA  A ,2) 
E RR=DRAW* SORES +34 R ERR 

ROC MON ,YR) =A (MON) f  BIMON) *P (MON , YR) +C (MON) *BRO+ERR 
3RO=RO (MON, YR) 

IF(RO(MON,YR)  .l_T.  j  .0)    RO  (MON,  YR)  =0  .  0 
I F ( MON  .  EQ  .  12 )    GO   TO  1013 
M  0N=M0M  +  1 
GO  TO  101*+ 

1013  CONTINUE 

C  COMPUTE   MEANS   AND   STAMOARO  DEVIATIONS 
C  OF  SIMULAT-D   VALUES  3Y    CALENDAR  MONTHS. 

DO   25ul  MON=l,12 

SUMP (MON) =  0 . 0 
"      SUMP2 ( 10N) =0.0 

SUMR(MON) =  0 • J 

SUMR2 ( iON) =0. G 

DO  2?uj  YR=i,NYRS 

3UMP(  MON)  =SUMP  (MOvi  ) +P  (MOM  ,  YR)   


250  1 


S  UMP2  {  iON)  =SUMP2  (1uN)+PM0N,YR)»P(M0N,YR)  ______ 

SUMR(MON) =  SUMRMON) +  R0M3N,  YR) 

SUMR2 ( HON) =SUMR2( HON) +  R0( MON , YR) *R0 (MO N , Y R ) 

5UMP2  (  IOiM  )  =  SCUT  (  (  5UMP2  (NGN  )  -  (  S  JMP(MON)  *SUMP  (HON  )  )  /  NY^S)  /  (NYR5-1)  ) 
5UMR2 ( M  ON  )  - 5QR  TJ  (  5  UHR2CM  ON)-(SUMR(M  uN)  *3GM  R  (  M  ON)  )  /NYRS)  /(NY  <  5-1 )  ) 

SUMP (MON) =SLMP{ NCN) /NYRS 


WRITL(6,lClo) 

r;Ll\  »  bUr  r  I  r  Jl\  J  f 

li*P2  (MC.N)  , 

{ YR,P ( MCN, 

IfPJ  , 

yR=UNVRSl 

1016 

FCRMATI*0  rai 
l.*,Fb. 4/(12(1 

N0^*,i3/* 

AVERAGE*  , 

F7«  3 

,*       STD.  DEV 

10 15 
1017 

MRU  E<  6*  1017  J 
FGkMAT(*0  rum 

•  *  j  i  J/ v  A 

I YR»RQ( MON 
V  Er-AGt*  »  F7 

» YR ) 

,YR=lt  NYRSI 
SFJ.  DEV.*, 

C  CCN 

lFo.^/i  12(  I  4 ,  P 

6.2))) 

AMS  OF  PR  EC  1  P  .  A 

o  Y     _  .Trio. 

DO  1C19  M0N=1 t12 

1 02  0 

DO    1021  YR=1, 

CLAS=CLAS* 
1  =  1 

1C23 

IF (P (MON  ,  YR  )  . 
1=1  +  1 

L  T  .CL AS  )   GO   TC  1 

C  __  ? 

CLAS=CLAS+CLASfc 

GO  TC  1023  _____ 

1C22 
102  1 

N INC  LI  I ,H0N3  = 
CONT  tNUE 

MNCL1I  ,MGN)+i 

1C19 


NT 


AT  0-00   INCH   AND  HAVE   A  CLASS 

i 

WIDTH  0F*,F5.2 

t  *   I NCH • * ) 

CG  1024  M0N=1, 

1C24 

WKi  TE ( 6  ,  1C25) 

M ON , ( I ,NlNCL(l»f*CN) » I = 1 

1025 

D GR A M~OF "RAINFALL  FOR  M 

□  NTH  N0-*Tf3/(10(I6,I4)  )  )~ 

DO   1G23    1  =  1 , 5C 
1028   N  I N  C  L  (  I  »  V  C  N  )  =  0 

Gt_  1C29  YR=1»NYRS 
CL  AS  =  CL  ASto 


IGjC    IF ( RG (MON, YR) .L T.CLAS)    GO  _JQ_  1 03JL  

1=1  +  1 

 CLAS=CLAS+CLASW 

GO   TO  1030 

1031  N1NCLCI  LHCN )  =  M  NC  L  ( I  ,  MON ) + 1 
1029  CONTINUE 

1027  CCNTINUE 

DO   1C32   MCN  =  1 ,12 

1032  WR1TE(6, 1^33)   MON,    (  I  , N  INCL ( I  ,M  JN  ) , 1  =  1,50) 

1033  FORMAT (*0  Hl^TCGRAM  OF   RUNCFF   F_r    WONT F   N0.*,I3/( 10(16,14))) 

IF(SET.Eg.NSETS)  GO   TO  1035     

StT=SET+l 

GO   TO   1C34   '  

1035   WRITE(6,1036)  RN,N1,N2 

1C36  FORK AT ( * C  LAST   RAN CON   NUMBERS  ARE*/3I17) 
IF  (STN.  Eg.STNS )   GO   TO  1043 

 5TN  =  STN  +  1  -  

GO   TC  10^4 
1043  STOP 
END 
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<;(  i  r:  r    ;i  |j  |  vc    r  a 

01 MF'iS I QN  TAB 

"pjtfofr  on" 
if( ient.ne .1  ) 

Kl=16777219 


JJ-IMTi  10.03!  Ft  QATIN2  )  /  <?  )  + 


on  1000 
00  loo  1 

tool  TABU,  J 
"1  rinn    r 0  • >  t t  r  •  •  1 


AR ( I  I , J  I  ) 


LL 


II  =  P  r  (  1 

Ml =M1 *K 1 


GO  T n  1004- 
0"    JJ=If  T  (  3  0  .  ( 

N2=N?*K1 
04  [  c  =  I c  + 1 

FNn 


SU BR OUT  IMF   R NM 

S  IM=0.0  ~ 

DC   tp   1  =  1,1? 

N0=N0*r6  77  7?l<> 
fU  =  FLOAT  (  NT  )  /  ? 


1 


S' 


Tt  IP  M 
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